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SUMMARY 


Analytical design techniques for active and passive control of aero- 
elastic systems are presented. These techniques are based on a rational 
approximation of the unsteady aerodynamic loads in the entire Laplace 
domain, which yields matrix equations of motion with constant coeffi- 
cients. 

Some existing rational approximation schemes are reviewed, the matrix 
Fade approximant is modified, and a new technique which yields a minimal 
number of augmented states for a desired accuracy is presented. 

The state-space aeroelastic model is used to design an active control 
system for simultaneous flutter suppression and gust alleviation. The 
design target is for a continuous controller which transfers some meas- 
urements taken on the vehicle to control command applied to a control 
surface. The control parameters are constant and they are optimized to 
minimize any desired combination of gust response parameters in a way 
that assures stability over the range of varying aerodynamic parameters 
in the entire flight envelope. 

Structural modifications are formulated in a way which enables the 
treatment of passive flutter suppression system with the same procedures 
by which active control systems are designed. 



INTRODUCTION 


Background 

In recent years, extensive research has been carried out to develop 
active control systems for controlling aeroelastic response. In these 
systems aerodynamic control surfaces are operated according to a control 
law which relates motion of the control surfaces to some measurements 
taken on the vehicle. The aerodynamic forces generated by the control 
surfaces modify the overall forces in some favorable way as defined by 
the performance requirements. 

The increasing emphasis on fuel efficiency and the advances in aero- 
dynamic and structural design techniques results in increasing payload to 
structural weight ratio. This increased structural efficiency results in 
lower elastic mode frequencies. Thus the elastic modes are more easily 
excited by atmospheric turbulence and pilot control inputs. Several 
development and flight test programs, such as the B-52 Load Alleviation 
and Mode Stabilization (LAMS) [1] and Drones for Aerodynamic and 
Structural Testing (DAST) [2], demonstrated aeroelastic control concepts 
for augmented rigid body stability, maneuver load control, ride control, 
fatigue reduction, gust alleviation and flutter suppression. The DAST 
program is still in progress and will continue for several years. The 
benefits of the aeroelastic control systems are increased aircraft 
maneuverablity, rider comfort, service life and flight envelope for a 
given structural layout. 
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Flutter suppression is fundamentally different from the other active 
control applications mentioned because the structural stability of the 
vehicle is Involved. Loss of a flutter suppression control system may 
result in almost immediate major structural failure and the loss of the 
aircraft. Consequently, the active control system hardware must be of 
high reliability, and structural stability should be assured for all 
possible structural configurations in the entire flight envelope, 
extended to include safety margins. Some recent analytical developments, 
wind-tunnel tests and flight-test demonstrations [2-6] show the potential 
and feasibility of active flutter suppression control systems. However, 
to apply modern control techniques efficiently for optimized systems, 
more refined aeroelastic modeling has to be carried out. 

The main difficulty in modeling an aeroelastic system for control 
analysis lies in the representation of the unsteady aerodynamic loads. 

The time lags in the development of these loads result in analytic ex- 
pressions which contain non-rational terms. Rational approximations are 
needed to obtain finite-order models which can be solved by the methods 
of linear algebra. 

Once a suitable approximation for the aerodynamic loads is chosen, a 
state-space matrix equation of motion can be constructed. Because 
aerodynamic loads are a function of flow conditions, each point in the 
flight envelope has a different equation of motion. A root-loci 
analysis yields the variation of the system frequencies of oscillation 
and damping ratios with the dynamic pressure for a given Mach number. 
After adding to the state-form equation, a control (input) term and a 
measurement (output) equation control, analysis is performed. 
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There is no all-purpose control technique. A suitable one is chosen 
in accordance with the performance requirements, available control 
means, precision of the mathematical model, measurement accuracy and 
on-board computer capacity. The control analysis given in this study is 
for a constant parameter, partial-feedback, continuous control system. 

A combination of modal and control surface response to atmospheric tui — 
bulence at various flight conditions is used as a cost function in such 
a way that stability is assured and a desired gust response parameter is 
minimized. 

Structural modifications are more traditional aeroelastic control 
means. Passive compensators such as mass balance or stiffness tuning, 
with the associated weight penalty, can be more practical than active 
means In many instances and should not be overlooked In the design 
process. Passive means are formulated in this work as input and output 
terms such that they can be included in the active control design. 

Survey of Literature 

The scope of this work goes across several major topics such as ae- 
roel astici ty , structural dynamics, turbulence, control analysis and op- 
timization. The literature survey given in this section is a selective 
review. Many other relevant books and papers appear in the literature. 

The idea of using an active control system to move aerodynamic con- 
trol surfaces for flutter suppression was pursued in the last decade by 
several aircraft manufacturers and research institutions. The B-52 CCV 
program was the first to demonstrate such a system in flight. This and 
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some other feasibility studies [7,8], wind tunnel tests [4,6] and flight 
tests [2,5,9] show the potential and feasibility of such systems. 

The structural analysis methods needed for earoelastic analysis are 
reasonably well developed. For example, the infinite dimensional spaces 
required to describe exact motion of aeroelastic systems can be reduced 
to finite dimensional spaces by the technique of truncated normal modes 
[10]. Further, normal mode analysis of the structure can be done by the 
well-developed finite element method [11]. 

The unsteady aerodynamic loads due to an oscillaning wind section in 
Incompressible, two-dimensional flow were first developed by Theodorsen 
[12]. Extensions of Theodorsen 's derivation are given by Timman and Van 
der Vooren [13] for subsonic flow and by Garrick and Rubinow [14] for 
supersonic flow. A thorough review of these and other derivations, in- 
cluding Possio's integral equation for the pressure distribution on a 
two-dimensional typical section due to small oscillations of assumed 
modes and the corresponding integral equation for three-dimensional 
wings, is given by Bisplinghof f , Ashley and Halfman [10]. Among the nu- 
merical techniques for solving the Integral equation for a three-dimen- 
sional wing in subsonic flow is the doublet-lattice method of Albano and 
Rodden [15]. 

The availability of analytical tools for calculating the aerodynamic 
influence coefficients for lifting surfaces in simple harmonic motion 
leads to the classical V-g method for predicting flutter boundaries 
[10,16,17]. The flight conditions at which the dumping g ceases to be 
negative defines the flutter condition. 
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Sears 118] developed convolution integrals for the unsteady 
aerodynamic loads due to arbitrary motion of a typical section in incom- 
pressible flow and generalized Theodorsen's function for arbitrary mo- 
tion. Edwards [19,20] applied the Laplace transformation to Sears' ex- 
pressions and generalized Possio's integral equation for arbitrary mot- 
ion. He also presented the aeroelastic equations of motion in a way 
that manifests the benefits of rational approximations of the aerody- 
namic loads> especially when an active control system is used. 

R. T. Jones 1 21] first introduced the use of rational Laplace trans- 
fer functions to approximate unsteady loads. The importance of such ap- 
proximations increased with the progress in active control of aero- 
elastic systems. Roger et ai* [3] used Pade approximants [22] separately 
for each term of the aerodynamic influence coefficient matrix. Roger [23] 
Increased the approximation efficiency by using common denominator roots. 
Vepa [24] and Edwards [19,20] used matrix Pad^ approximants to deal with 
all the Influence coefficients simultaneously. 

The analytical tools for calculating the response of flexible airplanes 
to atmospheric turbulence are summarized by Pratt [25]. The mechanism 
and theory of turbulence are given by Hinze [26]. Dryden [27] and von 
Karman [28] gave practical formulas for the statistical description of 
atmospheric gust velocity. 

Like any other component of a flight vehicle, an active control sys- 
tem must comply with the Federal Aviation Regulations [29]. An example 
of regulation considerations for an active flutter margin augmentation 
system is given by O'Connell and Messina [9]. 
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Aeroelastic active control design concepts are functions of the 
available aeroelastic modeling techniques. Nissim [30l developed the 
aerodynamic energy concept, uhich is based on the classical assumption 
of simple harmonic motion. This concept was applied by Nissim and 
Abel [4] and Abel a^. [31] . 

Rational approximations of the unsteady loads have been used 

to construct finite state-space mathematical models. Classical control 
techniques [31^321 and optimal control techniques [19,32J have been used to 
design flutter-suppression systems. The concepts and methodology of op- 
timal control theory are given by Bryson and Ho [341. Ashkenazi [35] 
applied some of these concepts to the constant parameter insensitive 
control of a system with variable dynamics. 

When control parameters are optimized with respect to a given cost 
function, and there is no direct algebraic optimal solution, a minimiza- 
tion procedure can be used, Fletcher and Reeves [36] introduce the 
conjugate gradients technique. Davldon [37] introduced a rapidly 
convergent descent method which has been modified by Fletcher and Powell 
[38] and by Steward [39] to include finite difference approximation of 
the gradient vector. 

The costs and benefits in using active control system to suppress 
flutter and alleviate gust response must demonstrate superiority over 
the more traditional passive flutter-suppression means. Structural op- 
timization procedures for satisfying flutter requirements are given by 
Markowitz [40] and Haftka and Starnes 1411. The decoupler pylon con- 
cept, which avoids violating the flutter requirements when an external 
store is added to the wing, was developed by Reed ejt al. [42]. 
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Report Outline 


In the Aeroelastlc Equations of Motion Section the typical section in 
incompressible and subsonic flow is introduced first. Unsteady aerodynamic 
loads in three-dimensional flow are discussed, and the use of truncated normal 
modes is presented. 

In the Flutter and Gust Response Section the V-g method and the root-loci 
approach for calculating flutter boundaries are presented. The continuous 
gust loads and the associated structural dynamic response are defined in 
statistical terms, and the basic relations are given. 

In the Finite State-Space Modeling Section the initial theoretical develop- 
ments original to this work are presented. Some existing techniques for 
rational approximation of the aerodynamic loads are reviewed, the matrix Fade 
approximant is modified, and a new minimum-state technique is developed. 
Numerical examples for approximating two-dimensional aerodynamics in incom- 
pressible and subsonic flow and three-dimensional aerodynamics in subsonic 
flow are given. 

In the Active Flutter Suppression and Gust Alleviation Section active 
control of aeroelastlc systems is treated. Design targets are defined, 
some existing techniques are reviewed, and a new approach for simultaneous 
flutter suppression and gust alleviation is developed and illustrated. 

In the Passive Flutter Suppression Section a new approach to passive 
means for flutter suppression is presented. The passive means are formulated 
as input and output terms of the active control equations. Numerical 
examples of stabilizing the typical section using a concentrated mass and 
a dynamic vibration absorber are given. 

In the final two sections. Concluding Remarks and Recommendations, the 
material presented in this paper is summarized and recommendations for 
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future research presented. 


Summary of Contributions 

1. Rational approximations of the unsteady aerodynamic influence 
coefficients are investigated, the matrix Fade technique is 
modified and a new "minimum-state" technique is developed. 

Numerical examples are used to demonstrate significant advan- 
tages of the minimum-state technique over other techniques. 

2. The minimum-state approximation is used to construct state- 

space aeroelastic control equations with flight-condition dependent 
coefficients. Statistical turbulence parameters are formulated 
as process noise. 

3. A control cost function is defined such that a desired gust 
response parameter is minimized and stability over the entire 
flight envelope is achieved, if at all possible with the 
available control means. 

4. A pole-assignment technique for a partial-feedback, zero-order 
compensator is modified to accommodate simultaneous pole as- 
signments at different flight conditions. 

5. Higher than zero-order compensators are augmented to the aero- 
elastic control equations, and the output equation includes an 
input term to accommodate acceleration measurements. 


6. Passive control means are formulated as Input and output terms 
such that they can be treated as a special case of active con- 
trol means. 
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SYMBOLS 


a 

ao 

®g ,i 

®u 

b 

bp 

c 

S 

Ct 

C(s') 

Do(s) 

f (xi ),g(xi ) 


g 

h 

Hn‘ (2) 
Im(*) 

J e 
J o 

J n(z ) 
k 

k c 
kp 


dimensionless elastic axis location (Fig. 1) 
speed of sound 

cost function Heights^ Eq. (5.16) 
control command participation in zj, Eq. (5.14) 
reference semichord 
viscous damping 

dimensionless control surface hinge location (Fig. 1) 
lift coefficient derivative with respect to angle of at- 
tack, La/2bq 

generalized Theodorsen function (Eq. (2.12)) 
open-loop characteristic polynomial 

coefficients of spacial longitudinal and vertical veloc- 
ity correlation, Eq. (3.11) 
structural damping 

plunge displacement (Fig. 1), altitude in Fig. 2 
Hankel function of order n 
imaginary part of • 

least-squares cost function, Eq. (4.2) 
global control cost function, Eq.(5.16) 

Bessel function of order n 
reduced frequency, «b/V 
integer defined in Eq. (5.5) 
spring stiffness 
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Kn(s') 
1 c 

L 

Ua 

Lg 

tn 


me 


m p 


ms 


M 

n 


He 


N 

P 

P" 

q 


Q i 
r 


ra^,rp^ 


Re(*) 

s 

s' 

t 

Ti 


modified Bessel function of order n 
compensator order 
lift per unit span 

lift per unit span due to unit pitch angle 
scale of turbulence, length units 

order of the denominator polynomial in a rational approx- 
imation 

number of output measurements 
concentrated mass 
uing mass per unit span 
Mach number 

pitching and hinge moments (Fig. 1) 

number of degrees of freedom 

order of control equation (5.1) 

order of Roger's approximation, Eq. (4.5) 

pressure difference across the wing 

ampl i tude of p 

dynamic pressure, ipV^ 

polynomial coefficients of Dq(s), Eq. (5.21) 
magnitude of s' = r exp(i0) 

radius of gyration of the typical section and the control 
surface, divided by b 
real part of • 

Laplace transform variable 
nondimensional ized Laplace variable, sb/V 
time 

constants defined in Eq. (A. 2) 
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u,U(s) 


u' , v' ,w' 


Ug,Ug 

V 

w" 


X1 »X1 ,Z1 


Xte 

Xa.Xp 

YnCz) 

z 

Z« 

ZdrZdCs) 


control command variable and its Laplace transform 
disturbance velocity components in Xi»vi and zi direc- 
tions, respectively 

vertical and longitudinal gust velocity 

airspeed 

amplitude of w' 

streamuise, spanuise and upwards spacial coordinates 
Xi coordinate of the trailing edge 
dimensionless center of mass offsets (Fig. 1) 

Bessel function of order n 

defined in Eq. (A. 4) 

z-i coordinate on the wing surface 

gust response design variable and its Laplace transform 


a 

Vj 

rcn) 

e 

S 

^ M e ( ^ ^ 

A 


P 

O'wg* 


angle of attack 

control surface actual and command rotations 
aerodynamic roots (Eq. (4.5)) 
gamma function 

error of least-squares approximation 
concentrated mass location. Fig. 18 
power spectral density of Wg 
disturbance velocity potential 
modal deflection 
wave length of sinusoidal gust 
mass ratio, ms/npb^ 
air density 

mean-square value of Ug 
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ffzd ,i^ 

6 

(j 


wh»wa» wp 


Wn 


n 


mean-square value of zj at point i of the flight envelope 
argument of s' 

vibration frequency» rad/sec 

plunge, pitch and flap uncoupled natural frequencies 
natural frequency 

aeroelastic eigenvalues defined in Eq. (3.2) 
modal damping 


Matrices 
ad j [ • I 

[adl, [Ool/ (avl 
[A(s') I 

[ Aap ] 

[AcI»fBcl 

lAj» ililBj, il 
( A 1 ] , I Bi ] 

[BI 
lBi ] 

IC], [Cel 
[C'l, [Cc'l 
[C(R',ki)l 
[ C 2 ] 

ID1,[D'1 
[ ] , [ O 2 1 
[El,[E'J 
[Eil 

lF(k)],[G(k)] 


adjoint matrix of I*] 

frequency response coefficients in Eq. (3.10) 
aerodynamic influence coefficients 
approximated value of lA(s')] 
compensator parameter matrices, Eq. (5.5) 
gust loads influence coefficients 
as defined in Eq. (4.12) 
known data matrices in Eq. (4.1) 
damping matrix 
defined in Eq. (4.15) 
control gains, Eq. (5.7) 
defined in Eq. (5.8) 
as defined in Eq. (4.27) 
defined in Eq. (5.12) 
defined in Eqs. (4.23) and (4.25) 
defined in Eqs. (4.15) and (4.18) 
defined in Eqs. (4.22) and (4.25) 
defined in Eqs. (4.15) and (4.18) 
real and imaginary parts of [A(ik)] 
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[FiJ dynamic matrix of aeroelastic control system^ Eq. (5.1) 

^g(s)^ as defined in Eq. (5.19) 

control input distribution vectors, Eq. (5.1) 
output matrices of Eq.(5.2) 

[Half (Hdl, [Hy] partitions of tHl, Eq. (5.2) 

[Hq] acceleration output matrix, Eq. (5.3) 

[ I ] unit matrix 

[K],[M] stiffness and mass matrices 

[KiIjIKzI defined in Eq. (4.15) and (4.22) 

{ La(t) La(s) aerodynamic load vector and its Laplace transform 
[Nil matrix polynomial coefficients in Eq. (5.20) 

[p^l coefficients of matrix polynomial approximation 

iQpt] piston theory limit, page 110 of Ref. 116] 

[R] approximation root matrix, Eqs. (4.9) and (4.25) 

[R'l defined in Eq. (4.24) 

^Ri ISi ] , [S 2 1 defined in Eq. (A. 2) 

[T] similarity transformation matrix, Eq. (4.16) 

^x^,^X(s)^ deflections in generalized coordinates and their Laplace 
transform 

^ x« ^ Xai ^ Xa 2 ^ aerodynamic augmented state vectors 
•{y }-,'{y(s) ^ measurement vector and its Laplace transform 

^z^,^Z(s)}- state vector and its Laplace transform 

compensator state vector 
modal deflection vector 
noise distribution vector, Eq. (5.10) 

Subscripts 

a aerodynamic 
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d 

design 

e 

related to least-squares problem, Eq, (4.1) 

eq 

equivalent 


1 values at which approximated and tabulated terms match 

i / j 

of i-th and j-th mode, respectively 

[• ,i»J 

row i of [•] 


column j of [•] 

• 

ij-th element of !•] 

1 

index of tabulated data 

mo 

maximum operation 

nc 

nonci rcul atory aerodynamic 

P 

passive 

r1,r2 

reference values of 1 in the least-squares equations 

s 

structural 

si 

sea level 

Superscripts 

t 

matrix transpose 

- 

average value 

Abbreviations 

C.G. 

center of gravity 

DOF 

degrees of freedom 

CCV 

control configured vehicle 

E.A. 

elastic axis 

FAA 

Federal Aviation Administration 

kPa 

kilo-Pascal, 10^ Newton/m^ 
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AEROELASTIC EQUATIONS OF MOTION 


Unsteady Aerodynamics 

The development of the partial differential equation for unsteady po- 
tential flou is presented in numerous textbooks and papers. An overall > 
thorough presentation is given by Bisplinghoff et ai. MOJ. The lineal — 
ized partial differential equation for unsteady* compressible flow reads 

1 ^ 2 $/ ^ 2 $/. 

^2^/ + 2V + V2 = 0 (2.1) 

ao^^St^ Sxi&t ^xi^. 

where $' is the disturbance velocity potential. The disturbance velocity 
components are 






I 

II 

> 

, w' = 

Sxi 

3yi 



where x^, yi and zi are the streamwise* spanwise and vertical coordi- 
nates, respectively, when flow over a wing is studied. The associated 
disturbance velocities u', v' and w' are assumed to be small compared to 
V. 

For the unsteady lifting part of the solution of (2.1), a thin wing 
is replaced by its "flat plate" projection on the Xi-y^ plane. The lin- 
earized condition of no flow across the wing provides the primary bound- 
ary condition of (2.1), 
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w'(xi,yi,0,t) 


(2,3) 


dZa dZa 

St SXi 

where Za is the local Zi coordinate of the -uing. 

The unsteady flow is antisymmetric with respect to the x^-yi plane. 
The pressure difference across the wing is 


p(xi ,yi , t) 


= -2p 


S s • 

V + — 

Sxi StJ 


Cxi , yi , 0* , t) 


(2.4) 


Kutta's hypothesis of smooth flow off a sharp subsonic trailing edge 
requires that p(X'te'yi't) = 0. Substituting xi = x-te in Eq. (2.4) gives 
another boundary condition for Eq. (2.1)» 


■ S 
V 


St. 


^'(Xie»Vl>0,t) 


= 0 


(2.5) 


In incompressible flow, Eq. (2.1) reduces to Laplace's equation* 

7^$' = 0 (2.6) 

The solution of Eq. (2.1) or (2.6) has traditionally been simplified 
by assuming that the wing is undergoing simple harmonic motion* thus re- 
moving the independent variable t. Further simplification can be made 
if two-dimensional flow is assumed, thus removing the independent 

variable y . 

1 

For oscillatory motion* w'(xi*yi *0,t) = w"(Xi *yi )exp(iwt) and 
p(xi*yi*t) = p"(xi ,yi )exp(i«t) . For three-dimensional flow* Eqs. 

(2.1)* (2.3) and (2.4) combine to the Integral equation* 
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n"(xi,yi;k) ir 

V 8ttJ, 




-K(xi-f,yi-7);k,M)d7)df 


(2.7) 


uhere f and tj are dummy variables corresponding to xi and yi , respec- 
tively, and k = wb/V Is the reduced frequency. 

Landahl [43] gives an expression for the kernel function for nonpla- 
nar surfaces in subsonic flow. To simplify the formulation, the 

discussion here is limited to planar surfaces, for which the expression 
for the kernel function is 


exp [ -i k (xi )/b 1 


(yi-7))2 


Mlyi-7)| 


exp[-ik|yi-T)|/bl 


W(Xi-f)Z + (1-n2)(y,-7))2 /T+u7^ 


fexp[-iku|yi-T)l/bl 


-iku|yi-T)l/b] T 

du 

(1+u2)3/2 J 


( 2 . 8 ) 


U1 


where Ui = 


rh/(xi-^)2+(1-M2)(yi-T))2-(xi-t) 
(1-m2) |vi-d| 


Knowing w"(xi,yi;k), Eq. (2.7) is solved for p"(xi ,yi ;k,M)/q. Al- 
though this equation contains singularities, a solution can be obtained 
by numerical methods such as the doublet-lattice method [15]. 

The Typical Section 

The typical section is a simplified aeroelastic system consisting of 
a rigid plate of unit span connected to the ground by springs. It has a 
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trai 1 ing-edge flap and is exposed to tuo-dimensional flow as shoun in 
Fig. 1. Assuming no structural damping, the system equation of motion 
is 

+ iKsJ^x^ = •{La(t)f (2.9) 

where the displacement vector •{x}'* = [h/b,a,B], the aerodynamic load 
vector ^La}-* = l-Lb,Ma,M(jl and iMg] and iKsl are the structural mass and 
stiffness matrices 


1 

Xa 

X(J 

xa 

ra^ 

ru^+X(} (c-a) 

X(J 

to^+Xf} (c-a) 



[Kc] = msb^ 


Oh^ 0 0 

0 ra^cjjjZ 0 

0 0 




The unsteady aerodynamic loads for incompressible, two-dimensional 
flow were formulated by Theodorsen [121. The solution of Eq. (2.6) 
yields the pressure distribution due to oscillatory motion in plunge h 
pitching and control surface rotation g. Integration over the chord 
gives down-load -L, pitching moment, and hinge moment Mg per unit span. 

Following Sears [181/ Edwards [ 19 ] showed that these results may be 
generalized for arbitrary motion. The Laplace transform of ^La(t)^ is 

^La(s)^ = q[A(s')]^X(s)f (2.10) 
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uhere s' is the non-dimensional ized Laplace variable s' = sb/V/ and the 
aerodynamic influence coefficient matrix is 


[A(s')l = 2b2 


lMncIs'2 + iBncJ+CCs') 




+ [Kncl+C(s' 


)^RiHSi3 

j 


( 2 . 11 ) 


In this case, C(s') is the generalized Theodorsen function 

Ki (s') 

C(s') = (2.12) 

Ko(s')+Ki (s') 

where K Is a modified Bessel function of order n. The function C(s') 
n 

is analytic throughout the s’ plane except for a branch point at the 
origin, which requires a branch cut along the negative real axis. The 
matrices of Eq. (2.11) and the series expansion of the Bessel functions 
are given in Appendix A. Equations (2.11) and (2.12) are identical to 
those given by Theodorsen [12] when s' is replaced by ik. 

The two-dimensional oscillatory aerodynamic loads in subsonic flow 
were formulated by Timman and Van der Vooren [13] and are tabulated in 
the Manual of Aeroel astici ty [431. Vepa [2k] and Edwards [20] general- 
ized Possio's integral equation (2.7) for arbitrary motion in two-dimen- 
sional subsonic flow. 


Finite Wings in Three-Dimensional Flow 

Because a flying vehicle is a continuous elastic structure, the exact 
solution of the equations of motion requires infinite dimensional spaces. 
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Finite element modeling techniques [11] assume that the structure can be 
described by a finite number of degrees of freedom (DOF). A finite-ele- 
ment model of a complex uing structure may result in hundreds of DOF, 
for which as many natural frequencies and their associated mode shapes 
and generalized masses can be calculated. However, for gl obal -structure 
aeroelastic analyses such as flutter and gust response, only a limited 
number of the lower frequency vibration modes is required [10]. 

The formulation of the preceding section is now generalized for finite 
wings in three-dimensional flow. The Laplace transform of the n DOF, 
open-loop, aeroelastic system equations of motion (for stability analysis) 
is 

j^[Msls2 + iBsJs + [Kslj^X(s)^ = q[A(s')]^X(s)[ (2.13) 

Modal coordinates are used, which means that a discrete displacement 
is assumed to be a linear combination of a set of assumed modes, 

n 

xd(xi,yi,t) = X xd ,i (tl^fi (xi ,yi ) (2.14) 

i = 1 

Mhen the assumed modes are dimensionless normal modes with natural 
frequencies Wn, i, and when modal -damping coupling is negligibly small, 
the structural matrices of Eq. (2.13) are diagonal and their elements 
are 


and 


Ms.ii = //ms(xi ,yi )>Pi2(x, 
S 

Bs.ii - iM» ,i i(Jn ,i » 
Ks ,i i - Ms ,i i^n ,i ^ 


(2.15) 
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As stated previously, there are well-established numerical me- 
thods for solving Possio's integral equation (2.7) for oscillatory mot- 
ion (s' = ik). Substituting w' = Uj"exp(iwt) and z« = b4'jexp(icjt) into 
Eq. (2.3) yields the non-dimensional ized oscillatory downwash due to 
motion (of amplitude b) in the j-th mode 

wj"(xi,yi;k) 

= ik^fj + b (2.16) 

V <>xi 

Solving Eq. (2.7) with the downwash Eq. (2.16) yields the j-th 
non-dimensional ized pressure mode pj"(xi »Vi ;k,M)/q. The aerodynamic in- 
fluence coefficient matrix for a given M is calculated by 

1 PPp j"(xi ,yi ;k,M) 

A,ij(ik) =- (xi ry-i )dxidyi (2.17) 

bJJ q 

S 

Extrapolation of the aerodynamic influence coefficients to the entire 
s' plane will be discussed later in this paper. 
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FLUTTER AND GUST RESPONSE ANALYSES 

V-g Method for Flutter Analysis 

The flutter boundary is the group of points of the flight envelope 
(in terms of airspeed and air density or altitude) at which the aero- 
elastic system is neutrally stable or, in other words, a disturbance re- 
sults in simple harmonic oscillations. Thus, even though the computa- 
tional tools were (and to a large extent still are) limited to oscilla- 
tory air loads, the flutter condition can still be found. The most 
popular technique for calculating the flutter conditions is the V-g me- 
thod, originated by Smilg and Wasserman 117). 

Because the effect of the structural damping terms in Eq. (2.13) is 
usually small, because it is difficult to calculate damping accurately, and 
because omitting damping is usually conservative, structural damping is 
neglected. For a simple harmonic motion of frequence m, Eq. (2.13) can be 
written as 

f 1 pb2 1 

iMsI [Ksl + [A(ik)lHx(ik)^ = •{ 0 ^ (3.1) 

I 2kZ J 

A non-trivial solution of Eq. (3.1) with non-zero velocity and real w 
is possible on the flutter boundary only. To find these conditions, an 
artificial structural damping g is introduced by replacing 1/w^ in Eq. 
(3.1) by 

n = (1+ig)/«2 (3.2) 
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Equation (3.1) becomes 


[KsJ 


pb^ 

si + [A(ik)I 

2k2 


^X(ik)^ = n{x(ik)^ 


(3.3) 


which is an eigenvalue problem to be solved for n eigenvalues 17;. Fre- 
quency of oscillations, air velocity and damping parameter are determined 
from each eigenvalue by 


w;=)/l/Re(fi;) , V;=w;b/k , g ; =Im(n ; )/Re (J1 ; ) (3.4) 

Given p and M, Eq. (3.3) is solved for several values of k. Curve 
fittings of the variations of g and w with V yield the flutter velocity 
Vf and frequency Wf at the point of g = 0, or when g equals some es- 
timate of the damping present in the actual structure. The equivalent 
flutter speed Vfeq is calculated by 

Vfeq = ^/(p/Psl) Vf (3.5) 

The resulting Vfeq(M,p) does not necessarily match the actual equiva- 
lent velocity, V©q = vp/pTi Mao. To determine the actual flutter point 
for a given M, the flutter calculations are repeated iteratively with 
varying p until Vfeq(M,p) matches Veq(M,p). The group of actual flutter 
points, for various values of M, combine to form the flutter boundary. 

An example of a flutter boundary and a practical flight envelope is 
given in Fig. 2. Three envelopes are plotted: the operating flight en- 

velope which the airplane should not exceed in normal operation, (veloc- 
ities up to Vmo and Mach numbers up to ]%o) . the design envelope the air- 
plane is designed for (maxium velocity Vd and maximum Mach numbers Md) , 
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and the safety margin envelope which Is defined by aviation regulations 

(1.2 V, and 1.2 M,). On the Figure the h and h, refers to the maxi- 
d d mo d 

mum operational and design altitudes. Federal Aviation Administration 
(FAA) certification of a commercial transport airplane, for example, 
paragraph 25.629 of Ref. [29], requires that the unfailed airplane be 
designed to be flutter free to 1.2 times the design speed (V^j/M^) and 
demonstrated by flight test to be flutter free up to 

The graph of Vf should be out of the safety margin envelope for any 
weight loading and external store configuration. If not, as in the 
example of Fig. 2, the airplane should be placarded for lower ^mo/Mmo 
or modified structurally. Another alternative is to operate an active 
flutter suppression system which will drive the flutter boundary out of 
the safety margin envelope. 
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EQUIVALENT AIRSPEED 



M n 1.2 M, 

mo d d 


MACH NUMBER^ M 


Figure 2: Example of flight envelopes and flutter boundary 



Root-Loci of the Aeroelastic Modes 


The progress In developing analytical expressions and rational ap- 
proximations for the aerodynamic influence coefficients for arbitrary 
motion now enable the use of root-loci analysis to determine the flutter 
boundary. The problem is to find the roots si which yield a non-trivial 
solution for Eq. (2.13), namely 

IlMsls^z + [Bsls? + [Ksl - q[A(s'Ol| = 0 (3.6) 

The main difficulty in solving Eq. (3.6) is in calculating [A(s')l. 

The solutions for [A(s')I usually contain non-rational 
terms [191 and require lengthy calculations. Iterative procedures 
[19,20] are needed to solve Eq. (3.6) for s?. Given M, Eq. (3.6) is 
solved for various q, and the root loci are plotted in the Laplace s 
plane. For stability, all roots should satisfy Re(si) < 0. The flutter 
dynamic pressure qf is the q for which one of the roots is purely imagi- 
nary. The flutter frequency is 

Wf = Im(si) when Re(si) = 0 (3.7) 

In order to solve Eq. (3.6) by the methods of linear algebra, [A(s')l 
has to be approximated by rational functions of s'. Rational approxima- 
tions will be discussed later in this paper. 

The root-loci method provides "correct" vibration frequencies and damp- 
ing ratios at any flight condition, while the V-g method is correct on 
the flutter boundary only. Moreover, because p is not predetermined, the 
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root-loci -flutter dynamic pressure represents an actual flight condi- 
tion. 

In design for active control of an aeroelastic system, the root-loci 
method is more advantageous because it fits modern control techniques, es- 
pecially when rational approximations of the aerodynamic loads are used 
to construct finite order and constant coefficient system equations of 
motion, as will be discussed later. 

Continuous Gust Response 

An airplane traveling through the atmosphere encounters turbulence, 
which imposes gust loads on the structure. The turbulence and the asso- 
ciated gust response are defined in statistical terms, usually with the 
assumption that the process is stationary in time. A rigorous descrip- 
tion of gust response analysis of a flexible airplane is given by Pratt 
[ 25 ]. A summary of the work on mathematical modeling of turbulence is 
given by Hinze [ 26 ] . 

The atmospheric gust waves can be described as a superposition of 
waves of all wavelengths. The basic assumptions regarding airplane-gust 
interaction are that the turbulence pattern does not change appreciably 
during the time the airplane passes through (known as Taylor's hypothe- 
sis) and that the gust field is uniform spanwise. The latter assumption 
is correct for most airplanes, except very large ones. 

Consider an airplane flying horizontally with velocity V into a sinu- 
soidal vertical gust of velocity amplitude Wg" and wavelength The 
gust downwash at a longitudinal station Xi is 
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exp(-ikxi/b) expCiwt) 


(3.8) 


Mg = Mg" 

where « = 2 tiV/A and 


k = wb/V 


Substituting Wg"exp(-ikxi/b) for w" of Eq. (2.7) gives an integral 
equation for the nondimensional ized gust pressure difference pg'Vq. 
Substituting Pg'Vq for Pj'Vq in Eq. (2; 17) gives the gust influence coeffi- 
cients Agj. Adding the loads due to a sinusoidal vertical gust veloc- 
ity of unit amplitude to Eq. (2.13) gives 


-(j2lMs] + i«[Bs] + lKs]-q[A(ik) 1 


^X(iw)}^ = (q/V)-{Ag(ik) ^ 


(3.9) 


which is solved for the frequency response in generalized coordinates 
X(iw). Any gust-induced structural frequency response Zd(iw) can be de- 
scribed as a linear combination of the displacement, velocity and accel- 
eration frequency response in generalized coordinates. 


Zd( i«) 


ladl+i«[av]"«^{aoJ 


{X(iw) }• 


(3. 10) 


Before calculating the statistics of the gust response, the statis- 
tics of the atmospheric turbulence need to be known. The turbulence 
equations are very complicated, and closed-form solutions are available 
only in some specific cases. Isotropy and large Reynolds numbers are 
usually assumed for atmospheric turbulence. However, even with these 
assumptions the theoretical solutions are limited, and the mathematical 
models are based mainly on empirical and interpolation formulas. 
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The largest turbulence eddies (Iom wave numbers) are of permanent 
character. They contain relatively low energy per unit mass» but the 
energy dissipation by viscous effects is also very small. The larger 
eddies produce smaller ones through inertial interaction, thereby trans- 
ferring energy to them such that the energy containment is increased 
with decreasing eddy sizes. At the same time, the viscosity effects be- 
come more and more important and, starting at some wave number, the en- 
ergy decreases through dissipation. 

The probability distribution of the gust velocity is assumed to be 
Gaussian, which implies that its statistical properties are defined by 
its power spectral density (PSD) The PSD is a function of the 

coefficients of spatial longitudinal and vertical velocity correlation, 
fCx'i) and g(xi), which are defined as 


Ug(0)Ug(X'|) Wg(0)Ug(Xl) 

f(xi) = , g(xi) = (3.11) 

Isotropy implies Oug^ = <^wg^ which defines the relations between the 
statistical properties of the longitudinal and vertical gust velocities. 
Based on experimental results, Dryden [27] suggested that f(xi) can be 
approximated, for xi I 0, by 

f(xi) = exp(-xi/Lg) (3.12) 

where Lg is a turbulence scale and depends on atmospheric conditions. 

As derived on page 174 of Ref. 1261, the resulting PSD of the vertical 
gust velocity is 
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a„g2 Lg 1+3(wLg/V)2 


(3. 13) 


IT V [ 1+(wLg/V)2]2 

For high values of a, §hs((<’) of Eq. (3.13) becomes proportional to 
(wLg/'V)"*. Von Karman l28 1 indicated that it should be proportional to 
(wLg/V)*®''® (knoun as the Kalmogoroft spectrum law) and gave an interpo- 
lation formula for the energy spectrum function which reads 

Owg^ Lg 1+(8/3) (O/We)^ 

§„g(«) = (3,14) 

H V [ 1 + (u/'We) * 1 ’ ’ 

r(5/6) V V 

where We = Vtt — = 0.747 — 

r(1/3) Lg Lg 

The PSD of the gust response is related to that of the gust velocity 
by 


$zcl(cj) = I Zd( iw) I ^^wgCw) (3.15) 

The mean-square of the gust response is 
00 00 

= /|ZdCi«) I ^§wg(«)d<<; (3.16) 

0 0 

A numerical example of gust response of the typical section of Fig. 1 
in incompressible flow with V=275 b/sec was performed. (This is the 
first place in this work where new material appears.) The semichord b 
is use here and in later numerical examples as a length unit. The 
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frequency response was calculated using Eq. (3.9) with the structural 
parameters of Table 1, the gust loads of Eq. (A. 11), and Dryden's PSD 
of Eq. (3.13), with a^g = 1. (b/sec)^ and Lg = 50b. The resulting 
pitch angle frequency response and the pitch angle PSD are shown in 
Fig. 3. 


TABLE 1 

Three DOF typical section structural parameters 


«h = 

50 

rad/seo 

a = -0.4 

xa = 

0.2 

ua = 

100 

rad/sec 

b = 1. 

X(l = 

-0.025 

W(j = 

300 

rad/sec 

c = 0.6 

II 

0.25 

IJL = 

40 



= 

0.00625 
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sec 


Cl 


w, rad/ssc 

b) Pitch Angle Frequency Response 


w, rad/sec 

a) Input Gust Spectrum. 



c) Pitch Angle Spectrum. 




Figure 3: Dryden's gust PSD and gust response of a typical section in 

incompressible flou. 



FINITE STATE-SPACE AEROELASTIC MODELING 


As discussed previously, in order to solve aeroelastic equations 
by the methods of linear algebra# the aerodynamic Influence coefficients 
have to be approximated by rational functions of s'. Such approxima- 
tions enable finite, constant-coefficient, state-space modeling which 
best fits the formulation of modern control methods. 

The form of an approximating function depends on the ranges of inter- 
est, required accuracy, exact data availability and approxima- 
tion efficiency in later applications. Flutter reduced frequencies are 
usually within k = 0.1 to 0.5 (intermediate range). In stability and 
gust response analyses, however, lower reduced frequencies (k-»0) are 
also important for the following reasons: (a) static instability (div- 

ergence) might occur, (b) the power spectral density of atmospheric 
gust velocities is much higher in the low frequency range, as shown in 
Fig. 3, and (c) rigid body modes may be included. When transient re- 
sponse is of concern, higher values of k (or s') are required. This 
range is beyond the scope of this study. 

As discussed in the preceding section, three-dimensional aerodynamic 
routines are well developed for oscillatory motion only (along the 
imaginary axis of the s plane). Rational approximations are designed 
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to fit these data and are then extrapolated to the entire s plane by 
replacing ik by s’. The adequacy of this extrapolation will be 
demonstrated later in this section. The higher the order of the 
approximation (the number of denominator roots) , the more accurate is 
the approximation, but the resulting state-space model is also of 
higher order (less efficiency) . The numerical examples in this and 
following sections consist of the aeroelastic systems given in Table 2, 
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TABLE 2 


Aeroelastic systems for numerical examples 



Typi cal section 

Typical section 

Research uing 

Figure 

1 

1 

4 

Mach number 

0. 

0.7 

0.9 

DOF 

3: hta,B 

3: h,a,B 

4: three wing modes 
and one control 
surface 

deflection mode 

Structural 

properties 

1 

j 

1 Table 
! Eq. (2.9) 

a=-0.5 , b=1. 
c=0.6 

Abel [6] 


Aerodynamic 
data sources 

Eq. (2.11), 
Appendix A 

Manual of aero- 
elasticity [44] 

Doublet lattice 
method [15], 
supplied by Irving 
Abel of NASA-Langley 
and normalized for 
unit generalized 
mass 

Tabulated k 

0., 0.1, 0. 15, 

0., 0.06, 0.1, 

0., 0.1, 0.3, 0.5, 

values 

0.25, 0.3, 

0.16, 0.22 

0.7, 0.9, 1.3 


0.5, 1.0, 2.0 

0.26, 0.3, 0.4, 




0.5, 0.64, 0.8, 




1.0 




























Figure 4; Research wing geometry (taken from Ref. [6D. 
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Linear Least-Squares Solution 


A basic tool in rational approximations is the least-squares techni- 
que [45]. Given le data points, the problem is to find [Xel which ap- 
proximates the solution of the linear system 

lAi] [Xel - [Bil 1=1, le C4.1) 

neXme meXPe HeXPe 

where [Aj] and [Bil are functions of the tabulated data at each point. 

A "best fit" is obtained by minimizing the cost function 

1 e P e 

Je = Z 1 1 ei C4.2) 

1=1 j=1 i=1 

where £l^jj = [Ai^il^Xe^j^ - 

The least-squares criterion is satisfied by the solution of the normal 
equations of the least-squares problem: 

’ 1 e i e 

Z [Ai]*[AiJ [Xe] = I [Ail^tBi] (4.3) 

. 1=1 J 1=1 

When applied to approximating aerodynamic data, an index 1 refers to 
a known reduced frequency ki. The desired accuracy in different fre- 
quency regions can be controlled by (a) spacing the data points, (b) 
weighting the data by multiplying lAi] and iBil by different parameters 
for different 1 indices, and (c) requiring the approximation to match 
the data at some points. Matching can be done by either constraining 
some of the free parameters (as will be shown later) or applying large 
weights to the data of these points. 
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The matrix Fade method and the minimum- state method lead to 
least-squares problems which involve different terms of the data 
matrices simultaneously. The modal coordinates on which the aero- 
d3naamic data are based should be normalized such that the different 
terms are of comparable order of magnitude. 

Term-by-Term Pad^ Approximants 

R. T. Jones 1211 first introduced the use of rational Laplace trans- 
fer functions to approximate incompressible unsteady loads. Based on an 
exponential approximation of Wagner's indicial loading function (loading 
due to a step change in angle of attack), Jones' approximation of the 
generalized Theodorsen function reads 

0.0075 0.10055 

C(s') 0.5 + + (4.4) 

s'+0.0455 s'+0.3 

The branch cut along the negative real axis required for the exact 
function, Eq. (2.12), is replaced by two poles there. Edwards 118] used 
Jones' approximation to model the typical section with a state-space 
model of order 8 (six structural and two augmented states). 

In a more general aeroelastic system, unlike the incompressible, two- 
dimensional aerodynamics (Eq, (2.11)) case, there is no longer a single 
non-rational function which can be factored out of the influence 
coefficient matrices. This complication led to a Fade approximant 
technique [3,24], in which each term of the influence coefficient matrix 
is approximated by a different ration of polynomials in s'. 
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In a Fade approximant of order m> a least-squares technique is used 
to find the numerator and denominator polynomials which best fit tabu- 
lated influence coefficients for oscillatory motion. Since each approx- 
imation root adds one state to the equation of motion, the dimension of 
an n DOF model becomes 2n+mn*. In a realistic problem, the contrib- 
utions of some of these roots are negligible and are omitted. In the 
design problem of Ref. [31. for instance, a system of 27 structural 
modes was analyed with a mathematical model of the order of 200. 


Roger's Approximation 

Roger [23] realized that the aerodynamic influence coefficients may 
be approximated more efficiently by using common denominator roots. His 
approximation is 


N [Pjis' 

[Aap] = [PoJ + [Pi ]s'+[P2ls'2 + I (4.5) 

nxn j = 3 (s'+7'j-2) 


where the values of Yj-z are selected to be in the reduced frequency 
range of interest. The real coefficient matrices [Pj] are found by set- 
ting s' = ik and using a least-squares technique for a term-by-term fit- 
ting of tabulated oscillatory influence coefficient matrices [A(iki)l = 
[F(ki)l + ilG(ki)] for various 1 indices. The parameters of the normal 
least-squares equations (4.3) for the ij-th terms are 


[All = 


1 0 -ki2 


ki 0 


ki^ 

k 

kiri 

k 1^+71^ 


ki^ 

kl^+7N-2^| 
kirN-2 


ki^+rN-z 


2J 
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(4.6) 


,i j 1^1 


fPo .ij-l 

[BiJ = 

and (X«l = ■ 

• 

l-G , i j (k 1 )■' 


• 



•Pn ,13- 


An augmented state vector is defined by 


s 

^X.i(s)^ = -{X(s)^ 

nxl s+(V/b)Vi-2 


(4.7) 


Equation (2.13) can nou be re-transformed to a state-space matrix equa- 
tion of motion for stability analysis 


■ " 


- 



• 


• ■ 

X 


0 

[I] 

0 

• • • 0 


X 



- 1 

- 1 

- 1 

- 1 


• 

X 


-[Ml [K] 

-[Ml IB] 

[Ml [P3I 

• • • [Ml [PnI 


X 

X a 3 

- = 

0 

III 

-(V/b)YiIIl 

0 


Xa 3 

• 


• 

• 


• 


• 

• 


• 

• 


• 


• 

. • 


• 

# 


• 


• 

XaN 


0 

m 

0 

-(V/b)YN-2lI 1 


X aw 

. • 


• 



- 




where 


and 


[Ml = [Msl-ipb^lPzl 
[B] = [Bsl-ipbVlPi 1 
[K1 = [Ks]-ipV2[Pol 


The state-space model is of the order of Nn. Most applications (such 
as Refs. [6] and 1461) used N = 6. Roger's method was applied to ap- 
proximate the aerodynamic data of the research wing of Table 2. A best 
fit with N = 3, 4, 5 and 6 was obtained with ri = 0.2, Ya = 0.4, Ys = 
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0.6 and = 0.8. The results for uniform weighting are typified by the 
curve fittings of and ^(±k) , which are shown in Fig. 5. Two 

other cases were analyzed* one in which the approximation is required to 
match the data at k = 0.3 and one in which the matching requirement is 
for both k = 0 and k = 0.3. The matching requirement is imposed by as- 
signing large weights to the relevant data points. The total squared 
errors for the different cases is given in Table 3. 


TABLE 3 

Roger's approximation errors and aerodynamic roots 

J 

Sum of squared errors, m^ 


N 

Un i form 

Match 

Match k=0. 

Aerodynami c 


weighting 

k=0.3 

and 0.3 

roots (-yi) 

3 

.3406 

.4463 

2.5664 

-.2 

A 

.1128 

.1485 

. 1650 

-.2, -.4 

5 

. 0270 

.0437 

. 0483 

-.2, -.4, -.6 

6 

.0065 

.0079 

. 0087 

“•2^ • A f “*.6# ’“•S 
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Im(A ) ,m Ira (A 




Figure 5; 


Roger's approximation of 
at M = 0.9. 


A^ ^ and A of the research wing 
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Matrix Pad^ Approximations 


The matrix Pade approximant technique was introduced by Vepa 12A], 
uas modified by Eduards [20] and is further modified here. The approxi- 
mant is 


|-1i 


[Aap] 


sMil-lR] 


iPils'^+lPzls'+lPal 


(4.9) 


The data used to determine the approximant consist of tabulated 
steady-state influence coefficients [F(0)1 and oscillatory influence 
coefficients [A(iki)] = [F(ki)l + i[G(k])] for various 1 indices. A re- 
quirement of matching the steady-state aerodynamics yields 

[P 3 ] = -[RHF(O)] (4.10) 

[Pi] and [ Pz ] are related to [R] by 

[P 1 I = [R]([F(kfi)]-[F(0) J)/kfi2+[G(kfi)]/kfi (4.11a) 

and [Pz] = [F(kfz) ]-[R][G(kfz) 1/kfz (4.11b) 

uhere kfi and kfz are selected from the tabulated k values. Equations 
(4.11a) and (4.11b) are only approximations at other ki values, for 
uhich 


[Ai .iJMr]* [Bi ,i]t (4. 12a) 

and [Az ,a]*[Rl* [Bz (4.12b) 

uhere [A1.1] = ( [ F (k 1) ]-[ F (0 ) 1 )/k - ( [ F (k ri ) I- [ F ( 0 ) ] )/k ri ^ 
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[A 2 .ll = [G(ki)l/ki - [G(kr 2 >]/kr 2 


[Bi .il = [G(kri)I/kM - [GCki)]/ki 
and [Bz ,il = [FCki)l - [FCkrzJI 

where kri and kpz are selected from the tabulated ki values. 

Edwards [20] used kfi = kfz = kf> which causes the approximant to be 
exact at kf. He chose kri -* <» which yields C [ F (k n ) ]- [ F(0) ])/k ri ^ -► 0 
and [G(kpi)]/kri [Qpfl- He also used Eq. (4.12a) only to construct a 
matrix least-squares problem to be solved for [R]'^. 

The matrix Fade approximation is improved when both Eqs. (4.12a) and 
(4.12b) are used for least-squares fitting, which yields Eq. (4.3) with 



[Ai, il* 


iBi, il^l 

[All = 

[A2,i1^ 

, [Bil = 

1 


The matrix Fade technique was applied to the research wing of Table 
2 (with the wing modes only) for different kf and kr values. The sum of 
the squared errors and the aerodynamic roots (eigenvalues of [Rl) for 
the different cases are given in Table 4. 

Best overall accuracy is obtained by taking kfi = kri = max(ki) and 
kfz - kpz = ki at which accuracy is most important (case 1 of Table 4). 
The comparison with other methods is shown in Fig. 8. To yield the ap- 
proximation to match the data at k = 0.3, one can set kfi = kf 2 = 0.3 
(rather than = kpi), but there is a penalty in deteriorating accuracy at 
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Matrix 

Fade 

TABLE 4 

approximation errors and aerodynamic roots 

Case 

kfi 

kfz 

k ri 

krZ 

Sum of squared 
errors, m^ 

Aerodynamic roots 

1 

1.8 

0.3 

1.8 

m 

.183 

-.136, -.946, -1.928 

2 

0.3 

0.3 

1.8 

1.8 

.703 

-.138,-1.584, -2.474 

3 

0.3 

0.3 

1.8 

■ 

.701 

-.136, -.946, -1.928 

4 

0.3 

0.3 

■ 

1 

. 974 

-.141, -.961, -1.632 


higher reduced frequencies (case 3 of Table 4). The matrix Fade 
technique yields a state-space model of order 3n (Ref. [19]), which 
reads 


X 

X 

Xa 


0 [I ] 

- 1 - 1 
-[Ms! [Kg] -[Mg] [Bg] 

(V/b)[P3]-[Pl'][Kg] [Pzl-lPi'HBg] 



(4. 14) 


where lPi'l = (b/V) [ Pi HMgl" ’ 


Minimum-State Method 

Roger's method and the matrix Fade method provide 
straightforward linear least- squares solutions. Inspection 
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of the associated state-space equations ('1. 8) and (4.14), houever, 
suggests that their efficiency, in terms of number of augmented states 
per given accuracy in approximating the aerodynamic influence coeffi- 
cients, can be improved. 

Another approach taken here involves analyzing the state-space equa- 
tions first. A general state-space representation of an n DOF aeroelas- 
tic system uith m aerodynamic augmented states for stability analysis 
reads : 


X 

* • 

X 

X al 

The target now is to define the matrices of Eq. (4.15) in a way that 
is most convenient for the aerodynamic approximation, and without reduc- 
ing the generality of Eq. (4.15). The structural states ^x^ and ^x^ 
represent the physical and measurable motion of the structure, in terms 
of which the structiiral proper t i es and the tabulated aerodynamic data are 
given. The aerodynamic augmented states will now be redefined to reduce 
the number of free parameters in Eq. (4.15), without changing the struc- 
tural states. 

The eigenvalues of iRil represent the aerodynamic time lags. When 
these eigenvalues are constrained to be real and negative with linearly 
independent eigenvectors, there exists a similarity transformation ma- 
trix [Tl (Ref. 1471) such that 


0 

[I] 

0 


X 

[Ki ] 

[Bi ] 

[Di ] 


X 

[Ell 

[Ej] 

[Ri ] 


Xai 


n 

n (4.15) 

m 
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[RI = [T]-’[Ri J[TJ 


(4. 16) 


where IR] is diagonal with the same eigenvalues as lRil. By defining a 
new augmented state vector, 


^Xa2 ^ = ITI" ’-{xai 
equation (4.15) becomes 


^ • 

■ 

X 


X 

. r 

Xoz 


• ¥ 

*- 


0 

[II 

0 

[Kil 

[Bi ] 

[D2I 

[Ea] 

[E4 I 

[R] 


T 



X 


X 


XaZ 
» ¥ 


(4. 17) 


(4. 18) 


where [DzJ = [Di][T] , IE3] = [TI'MEi] 

and [E^ 1 = [Tl'^ [Ez J. 


A new augmented state vector is now defined as 

]xa^ = -{xaz^ + iRj-^EaHx^ (4.19) 

Differentiating Eq. (4.19) and using the last row of Eq. (4.18) gives 


^Xaf = [Eal^x^ + 


[E^]+[R] IE3] 


^x^ + [R]{xazK 


(4.20) 


Substituting Eq. (4.19) for ■{xaz^ gives 


^Xa^ = 


- 1 


lE^]+(R] [E3] 


^x^ + (Rj^Xa^ 


(4.21) 


Using Eqs. (4.19) and (4.21) to eliminate ^Xaz}’ from Eq. (4.18) gives 
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» • 
« 

X 

m « 

X 

X. 


0 

[I] 

0 


■ 

x 

iKz] 

lBi] 

[Dj] 


X 

0 

IE] 

IR] 


Xa 


where [Kzl = [K-iI - [D 2 r’[E 3 l 

and IE] = [Eh 1 + IR]-’ [E 3 I 


(4.22) 


The matrices of the second row of Eq. (4.22) are now defined such 
that the structural properties are separated from the aerodynamic ones. 

This yields 


X 


• X 
Xa 


0 m 

- 1 - 1 

-iMs+Mal iKs+Kal “iHs+Ma] iBs+Bal 

0 IE] 


• 



0 


X 

- 1 


. 

[Ms+t 1 al [D] 

. 

X 

IR] 


Xa 

- 


. . 


(4.23) 


Comparing Eq. (4.23) with Eq. (2.13) shows that the aerodynamic in- 
fluence coefficient matrix is approximated by 


lAap(s')] = [Pils '2 + [pjis' + [P3] 

nxn 

"j-l 

+ ID'] s'II)-lR'] lE'Js' 
nxmL mxm mxn 


(4.24) 


where the aerodynamic parameters of Eq. (4.23) are related to those of 
Eq. (4.24) by 

In,] = -Spb^lPi] , iBa] = -ipbViPz] , 


IK,] = -ipV2lP3] , ID] = ipV^lD'] 


(4.25) 


IE] = IE'] , [R] = (V/’b)[R'] , s = s'V/b 


51 



The approximation of Eq. (4.24) is constrained to match the influence 
coefficients for k = 0 and k = kf. Thus, 

iPlJ = C[F(0) J-lF(kf)J)/kf*+lD'](kf2[II+lR']2)-’[E'I (4.26a) 

[Pj] = lG(kf)l/kf+lD'l(kf2lIl+[R'F)-’[R'HE'] (4.26b) 

and 

[PsJ = [F(0)] (4.26c) 

For the other tabulated reduced frequencies, Eqs. (4.26a) and (4.26b) 

become approximations. Thus one obtains 

[D'][C(R',ki) ][El«([F(ki) ]-[F(0) I)/ki2-(lF(kf) ]-[F(0) ])/kf2 (4.27a) 

and [D'HC(R',ki) J[R'HE'] [G(kf)]/kf - [G(ki)]/ki (4.27b) 

where I C(R' ,k i) ]=(k [ I 1+lR' l^)- '-(kfM I I+[R' 1^ ) ' ’ 

Equations (4.27a) and (4.27b) constitute a non-linear least-squares 
problem to be solved for [D'J, (R'J and [E'J. The real part of the ta- 
bulated data is weighted in Eq. (4.27a) by 1/ki^ and the imaginary part 
is weighted in Eq. (4.27b) by l/kj. These are only examples which il- 
lustrate weighting the data for higher accuracy at low reduced frequen- 
cies. With uniform weighting Eq. (4.27) becomes 

ki2[D'llC(R',ki) ][E'I*:[F(ki) ]-rF(0)I-([F(kf) ]-[F(0) ])(ki/kf)2 (4.28a) 
and ki(D'l[C(R',ki)][R'l[E'I«‘[G(kf)lki/kf-[G(ki)] (4.28b) 

The minimum-state procedure goes as follows: 
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1. Set an initial number of augmented states m. 

2. Set an initial diagonal [R'] with distinct negative elements. 

3. Set an initial [D'J with rank [D'J = min(m,n). 


4. Form a linear least-squares problem to be solved for tE'J. 

For the uniform weighting of Eq. (4.28)» the parameters of the 
normal equations (4.3) are 


iXeJ = [EM 


[All 


■ki2[DM[C(RMki) 1 
.kilO'JtC(RMki) llRM- 


and [Bil 


'lF(ki) 1-[F(0) l-(F(kf) 1-[F(0) ])(ki/kf)2- 
.[G(kf)]ki/kf-[G(ki)] 


(4.29) 


The summation is done over all the data points except those 
corresponding to k = 0 and k = kf. 


5. Form a linear least-squares problem to be solved for [D']. 
The parameters of Eq. (4.3) are: 


[XeJ = [DM^ 


[All 


■ki2[EM*[C(R',ki) ] 
-ki[E']*[R'][C(RMki]). 


and [Bil 


[F(ki) J*-[F(0) I*-([F(kf) I-lF(O) ])*(ki/kf)2- 
-[G(kf)]*ki/kf-[G(ki)]* 


(4.30) 


6. Calculate the least-squares performance index Je by Eq. 
(4.2). 
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7. Repeat steps 4 to 6 to convergence. 

8. Use a minimization procedure to modify [R'l. The procedure 
used in this work is based on a modification of Davidon's me- 
thod [37] as given by Stewart [39l* 

9. Repeat steps 4 through 8 until the performance index converges 
to a global minimum. 

10. Calculate the approximant of Eq. (4.24) for the tabulated re- 
duced frequencies and compare with the tabulated data. 

n. If the accuracy of the approximant is insufficient in a lim- 
ited frequency range only, try changing the weighting of the 
data. If the insufficient accuracy is related to a specific 
mode, try changing the modal normalization to increase the 
weight of the associated terms. Repeat steps 4 through 10. 

12, If the accuracy of the approximant is still not satisfactory, 
increase m and repeat steps 2 through 11. 


Numerical Examples Employing the Minimum-State Method 

The minimum-state procedure is first applied to the typical section 
of Table 2 in incompressible flow. The tabulated aerodynamic influenc 
coefficients were calculated for pure imaginary s' values, using Theo- 
dorsen's formulation as given in Appendix A. The matching reduced fre 
quency is kf = 0,25 and the least-squares weighting is that of Eq. 
(4.27). The results are typified by Ai^2 which is, divided by -2b^, 
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the unsteady lift-curve slope A comparison betueen obtained 

a 'a 

from Theodorsen, Jones' approximation (Eq. (4.41), and the minimum-state 
method of order 2 (Eq. (4.24)1, is given in Fig. 6 for various values of 
s' = r exp(i6). 

The optimized aerodynamic roots obtained by the minimum-state method 
are -0.04746 and -0.2285, while those of R. T. Jones are -0.0455 and 
-0.3. While both approximations yield state-space equations of the same 
order (8), the minimum-state procedure is seen to yield a better approx- 
imation in the entire range of interest lr| < 1. and 60® < 9 < 120® 

Both approximations deteriorate as the branch cut along the negative 
real axis is approached. 

Although good accuracy is shown in this example for a wide range off 
the imaginary axis, one should not conclude that a rational approxima- 
tion, based only on oscillatory aerodynamic data, will always have sim- 
ilar accuracy in the entire s plane. As indicated by H. Ashley and W.N. 
Boyd, in a paper presented at Colloquium Honoring Hans Georg Kussner (1980), 
serious inadequacies might occur when the approximations are applied to 
three-dimensional wings in compressible flow. 

The minimum- state procedure, with the least-squares weighting of Eq. 

(4.27) and kf = 0.26, was also applied to the typical section of Table 2 
in subsonic flow of M = 0.7. The tabulated aerodynamic matrices were 
taken from Ref. [44]. The sum of the squared errors and the aerodynamic 
roots for various approximation orders are given in Table 5. The ap- 
proximants of order 1, 2, 4, and 6 are compared to the tabulated oscilla- 
tory in Fig. 7. The first-order approximation is quite poor. 

Starting with m = 2, it is up to the designer to make the trade-off be- 
tween the approximation accuracy and the order of the resulting model. 
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H- < 



Note: The values of r at the +'s along the curves for 9 = 60° 
and 120® are the same as for 8 = 90®. 


Figure 6: Rational approximations of of a typical section at M=0 as 

a function of s' = r expCiG). 
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TABLE 5 

Minimum-state approximation errorSf typical section at M=0.7 


Approxi- 
mation 
order (m) 

Aerodynamic roots (diagonal of [R']) 

1 

Sum of 
squared 
errors, b** 

1 

-.0600 

1791.0 

2 

-.0301,-. 1248 

73.3 

3 

-.0286,-. 1214,-. 9389 

36.0 

4 

-.0281,-. 1214,-. 9389,-3. 4279 

27.3 

5 

-.0260,-. 0594,-. 1501, -.791 3,-1 .969 

9.7 

6 

-.0099,-.0483,-. 1424, -.2727, -1.0241, -1.51 07 

5.6 


The last example chosen was the research wing of Table 2. The mini- 
mum-state procedure, with uniform data weighting (Eq. (4.28)) and kf = 
0.3, gave satisfactory results with four aerodynamic roots (m = 4). The 
optimized aerodynamic roots are -0.181, -0.384, -1.496 and -2.033. The 
sum of the square errors is 0.0332. This error is comparable to the ei — 
rors obtained with Roger's approximation with N=5 (Table 3). However, 
while the minimum-state approximation with m=4 requires four augmented 
states in the aeroelastic model, Roger's approximation with N=5 requires 
12 augmented states. 
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t1i n imum-state approximation of Cj^ of a typical section at 
M=0.7. “ 
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The minimum-state error is much smaller than the matrix Pade errors 


of Table 4. Curve fittings for selected aerodynamic terms using Roger's 
approximation with N=6, the modified matrix Pade approximant (case 1 of 
Table 4) and the minimum-state approximation are shoun in Fig. 8. The 
minimum-state approximations are within 3% of the data in the low fre- 
quency range (k < 0.5) and 5 % in the higher range. 

The minimum-state approximations will be used in the numerical exam- 
ples discussed In the following section. 
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ACTIVE FLUTTER SUPPRESSION AND GUST ALLEVIATION 


Problem Definition and Design Strategy 

The open-loop state-space aeroelastic equations of motion developed 
previously form the basis for control analysis. The purpose of this 
section is to demonstrate the use of some modern control techniques for 
flutter suppression and gust alleviation. 

The structural properties and the variation of the aerodynamic influ- 
ence coefficents over the flight envelope are assumed to be accurately 
known* and the aeroelastic system is assumed to be accurately repre- 
sented by a finite-order state-space equation of motion. The system is 
subjected to random excitation by a vertical gust, which is defined in 
statistical terms by its mean-square velocity and power spectral density 
(PSD). The frequency dependent aerodynamic loads caused by the sinusoi- 
dal gust are known over the flight envelope. 

The control means are aerodynamic control surfaces whose commanded 
rotations serve as input parameters. An arbitrary number of sensors 
measure discrete displacements, velocities or accelerations on the lift- 
ing surface. The measurements are assumed to be perfect. Compensator 
transfer functions close the loop by relating the control commands to 
the measurement signals. This is a constant parameter control system 
with no flight condition measurements. 
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A scalar design cost function is defined as a ueighted sum of 
mean-square gust response parameters at selected points of the flight 
envelope. The compensator parameters are optimized to yield a minimal 
cost function. Careful choice of the cost function parameters leads to 
a control lau, if there is any> uhich stabilizes the system throughout 
the entire flight envelope. The considerations in designing the cost 
function will be discussed subsequently. 

Optimal control theory [34,45] supplies a straightforward optimal 
solution for a controllable and observable linear control system with 
fixed dynamics and quadratic cost function. The controller and observer 
dynamics are uncoupled (the separation theorem), and their separate so- 
lutions combine to a compensator of the order of the system. In our 
case, the dynamics vary over a wide range, and the separation theorem 
does not hold. Furthermore, a compensator of the order of the system is 
not of particular significance. The analysis is expected to yield a 
lower cost function as the compensator order is increased. On the other 
hand, increasing order means higher analysis, hardware and debugging 
costs and lower reliability. A practical approach is to design for the 
lowest order compensator which stabilizes the system over the entire 
flight envelope, meets the control system hardware limitations and 
yields satisfactory gust response. 


The Control Equations 

A block diagram of the closed-loop aeroelastic control system is 
given in Fig. 9. The open-loop equation of motion (4.23) is now 
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supplemented with a control input term. The formulation and later nu- 


merical examples are for a single-input control system. The extension 
to a multi-input system is straightforward, except for the pole assign- 
ment technique which is described later in this section. The plant 
equations with no process noise are 

+ •{Gi}-u , nc = 2n+m (5.1) 

ncX 1 



'X ' 


' 0 

where ]z[ = 

X 

lx a- 

> "I gU = 

[Ms+Mal-^Gz^ 
. 0 


[F,] = 


0 [I] 0 

- 1 - 1 - 1 
-[Ms+Mal [Ks+Kal -[Ms+Mal iBg+BaJ iMs+Mal IDJ 


[El 


[R] 


The state vector and the dynamic matrix [Fil were discussed 

previously. The control input u is the commanded value of the con- 
trol surface rotation (Be in Fig. 1). The actual control surface motion 
states are part of ^z}^. It is assumed that the actuator is a linear 
spring with no damping or delay, such that is a function of the ac- 

tuator stiffness only. More complicated dynamics of the actuator system 
can be modeled by additional states or compensated for in the final con- 
trol law [ 6 ], While lFi] may vary considerably with Mach number M and 
dynamic pressure q the only term in •{Gi}^ which is a function of M and 
q is [Mai, which is usually small with respect to [M 5 ]. 
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Sensors located at various points of the structure take measurements 
uhich combine to form the output (measurement) equation 


+ -{j^u /fflc i He (5.2) 

mexi ncxi 

uhere [H] can be partitioned as 

[H] = j^EHdh [HvJ» iHal 

The vector ^y}- consists of discrete structural motion measurements, 
uhich are linear combinations of the modes uhich serve as structural 
generalized coordinates. When the measurements are displacement or ve- 
locity related, [H] is independent of M and q, [HaJ = 0 and = 0. 

The most practical measurement devices, houever, are accelerometers for 
uhich the measurement vector becomes 

• « 

^y^ = (5.3) 

uhere IHq] is defined by the modal displacements at the measurement 
points. Equation (5.1) is used to describe ■{y}’ of Eq. (5.3) in the form 
of Eq. (5.2) uith 


[HI = lHol[Ms+riaJ [^-[Ks+Kah-lBs+Bal, [Dij (5.4) 

-1 

and = [HoHMs+Mal 

Nou both [H] and •{j[ are functions of f1 and q. 
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A compensator of the order of Ic adds Ic states to the system. Fol- 
louing the formulation of Ref. [35], the compensator states are de- 
scribed in a compensator canonical form as 


■|zc^ - [Acl'jzc^ + [Bcl'jy^ 
IcXl 1 cxmc 


(5.5) 


where 


[Ac] 


[Aci I 

0 

0 


0 0 


[ A c2 I ^ 

0 lAckcl 


with [Acil 


0 1 
321-1 ®2 i 


[Ackcl = aic for odd I c 


and 


1 c^2 for even 1 c 
(1 c+D/2 for odd 1 c 


Equations (5.1), (5.2) and (5.5) are combined to form the equations for 
the augmented system. 


» ■ 
z 


[Fi I 0 


• • 
z 


• ■ 
U1^ 


• = 




■ + ■ 


Zc 


[BcJ[H] [Acl 


Zc 

. ^ 


[BeJ^jf 

. 


When the loop is open (u=0), the plant and the compensator are obvi- 
ously uncoupled. The loop is closed by a control law of the form 

u = [C]'{y^ + [Ccl^Zo^ (5.7) 

where [Ccl — [[Cci],[Cc 2 ]****>[('ckc]] 
with [Ceil = [0 1 ] and [Cckd = 1 for odd Ic 
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An extension of [CcJ to a multi-input system is given in Ref. [34]. 


Subst i tuti ng Eq. (5.2) into Eq. (5.7) gives 
u = [C'][H]^2^ + [Cc'J^ZcK 


(5.8) 


where 


1 

[C'l = IC] and 


1 

[Cc'l = [Ccl 


Substituting Eq. (5.8) into Eq. (5.6) gives 


z 


z 


c 


IFi]+^g4[CMIH] ^Gi^CcM 

[Bc] + lBcl^ JMc'HH] [Ac 1 +[Bc]{ JMCc' 


z 



(5.9) 


This is the closed-loop system equation of motion with no process 
noise. For no flutter in the flight envelope, all the eigenvalues of 
Eq. (5.9) should have negative real parts over the entire range of M and 
q. For a given system with me measurements, the free design parameters 
are the me gains in [Cl, compensator order 1 c» dynamics parame- 
ters ai through aic and the 1 exme parameters in [Bel. 

The cost function and the parameter optimization are discussed next, 
followed by discussion of the cases of zero-order and second-order 
compensators. 


Cost Function and the Optimization Procedure 

The cost function for which the free control parameters are optimized 
is based on the airplane response to atmospheric turbulence, which was 
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discussed earlier. In modern control terminology, the vertical 
component of the turbulence-induced gust velocity can be described as a 
process noise, as shown in the block diagram of Fig. 9. The plant Eq. 

1 5. 1 ) becomes 


where 


■|zf = lFi ]-{z^ + + -{rfwg 


f ° 

= {q/V[Ms+Mal'MAgCik) 

I 0 




C5. 10) 


As discussed previously, Wg Is defined in statistical terms by 
its mean-square a„g^, its power spectral density (PSD) §wg, and zero 
mean value. As shown in Fig. 3, the PSD is far from that of white 
noise. The distribution matrix in addition to being a function of 

M and q, is a function of the frequency w. The closed-loop gust re- 
sponse is calculated by the continuous gust response equations, modified 
to include the rational approximation of the aerodynamic influence 
coefficients and the control terms. 

The process noise term of Eq. (5.10) can be added to Eq. (5.9), which 
can be then solved for the gust frequency response of the whole system. 
Advantage is taken, however, of the special characteristics of the sys- 
tem, such that a much lowei — order system of equations is solved. Since 
the measurements are of structural motion only, vector ■{y}^ of Eq. (5.2) 
can be defined as 




(5.11) 
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To relate the frequency response of u to that of the variables 

of Eqs. (5.5)/ (5.7) and (5.11) are replaced by their frequency re- 
sponse, uhich yields 

U(iw) = lC2l^X(iw)^ (5.12) 




r 1 

• 

where 

[Czl = 

[Cl+[Co] iw[II-[AcI [Bel 

[ Hdl + i«I Hyl-w^ [ Ho 1 


From the second partition of Eq. (5.10), with the matrix definitions 
of Eq. (5.1), we obtain the equivalent of Eq. (3.9) for the closed-loop 
gust frequency response of the structural degrees of freedom. This 
equation is 

[-«2[Ms+Mal+tKs+Ka]+io(Bs+BaJ-iw(D] [iwll I-[R]1 [El 


-■{ G2 1 


^X(iw)^ = (q/V)^Ag(ik)f 


(5.13) 


Equation (5.13) consists of n equations to be simultaneously solved 
for ^X(iw)^. Calculating the frequency response from Eq. (5.10) would 
require solving 2n+m+mc equations simultaneously. The design frequency 
response of Eq. (3.10) is now modified to include a control term. 


Zd(iw) 


[adl + iwlayl-w^ I ao J 


■{x(iw)^ + auU(iw) 


(5. 14) 


where U(iw)is calculated from •{x(icj)}^ using Eq. (5.12). 
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The design mean-square frequency response is calculated by Eq. 

(3.16). Referring and Zd^i<^) to a point i in the flight envelope, 

defined by 11 and q, Eq. (3.16) becomes 

00 

Ozd ,i* = XI Zd ,i (iw) I ^♦wg(w)du (5.15) 

0 

The integration is done numerically, using the trapezoid rule [48], from 
zero to a frequency 505{ higher than the highest imaginary part of the 
eigenvalues of Eq. (5.9). 

The mean-square response of Eq. (5.15) serves as a local cost func- 
tion. A global cost function is defined as ueighted sum of local cost 
functions at selected points of the flight envelope 

Jo - Z ag^iCTzd.i^ (5.16) 

i 

The free control parameters, which were discussed in Section 5.2, are 
optimized for a minimal global cost function Jo by using Davidon's mini- 
mization method [ 37 ], as modified by Stewart [391 to accept different 
approximations of derivatives. 

The number of flight envelope points participating in Jo is limited 
by computer cost, as a gust response function evaluation is required in 
each step of the numerical integration of Eq. (5.15). The flight enve- 
lope points should include: 

1. Normal operation points (see Fig. 2) to minimize turbulence- 
induced vibrations and control activity during cruising 
f 1 ight . 
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2. Points within the safety margin envelope at which structural 
or control hardware limitations are exceeded. 

3. Points which might become unstable during the optimization 
process. 

It is clear that the points of the second and third categories above may 
not be anticipated in advance. After an optimization is performed for 
some cost function, root loci and critical response parameters are cal- 
culated for more points in the flight envelope. If the results are not 
satisfactory, the cost function is modified, and the optimization is re- 
peated, starting with the previous optimized control parameters. 

Gust response is infinitely large for points with pure imaginary 
roots. At each optimization cycle the initial control parameters should 
stabilize the system at all the cost function points. If the control 
means are sufficient, the optimization procedure tends to move roots of 
cost-function points, which are very close to the s-plane imaginary 
axis, away to the left. Thus, whenever a root-loci branch crosses to 
the right-hand plane within the flight envelope, a point close to the 
crossing point (but still stable) is added to the cost function. The 
weights ag assigned to the stability-related points, are relatively 
small to avoid unnecessary performance deterioration at normal operation 
points. 

When the control system is unable to provide satisfactory results, 
the design can be repeated with a higher order compensator. If the im- 
provement is not significant, a change in the number of sensors or their 
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locations^ or in the size and location of the control surfaces^ is 
needed. Sensor and control surface design> however, is beyond the scope 
of this study. 


Pole Assignment 


When a single-input, multi-output aeroelastic system is controlled by 
a direct, parti al -feedback control law (zero-order compensator), the 
only free parameters are the gains IC] of Eq. (5.7). The optimization 
procedure can now be applied to find the "best" control gains. Better 
insight and more efficient use of the optimization routine may be 
achieved, however, when assigned poles replace the control gains as 
free parameters in the optimization. 

The Laplace transform of Eqs. (5.1) and (5.2) yields the input-output 
transfer function 

^Y(s)}^ r T"1 

= [H] s[I]-lFi] (5.17) 

U(s) L J 


Substituting Eq. (5.7) with [Cel = 0 into Eq. (5.17) and premultiplying 
both sides by [Cl gives 


or 


where 


[Clj^[Hl[sm-[Filj = 1 

ICl [[H]-{g(s) ^+Dq(s)^ = Do(s) 


(5. 18) 


(5.19) 


^g(s)^ = 


ad j 


s[I ]-lFi ijJ^G^ » Do(s) = |s[IJ-[Fi]l 


In polynomial form 
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and 


■{g(s)^ = 


’Me nc- i 

I [Nils 
. 1=0 




ric nc- i 
Dots) = Y. tlis 
1=0 


(5.20) 


(5.21) 


An efficient algorithm for calculating the polynomial coefficients of 
•|g(s)^ and Dots) is given by Kailath [49]. The algorithm reads 


[Nil = [I] 



Qi = -(1/i)trace([Nj][Fi 1) 

for i = 1,nc 

(5.22) 

[Ni ] = [Ni-1 ][Fi ] + GU-i[I] 

for i = 2, no 



Once the polynomial coefficients are calculated for a point in the 
flight envelope, one equation (5.19) is constructed for each assigned 
closed-loop pole at this point. A complex pole assignment gives two 
equations, one for the real and one for the imaginary part of Eq. 

C5.19). A set of mo assigned poles gives mo equations to be solved for 
[C]. Different poles can be assigned for different flight conditions. 

Once the control gains are found, the closed-loop poles at any flight 
condition can be found by solving the eigenvalue problem 


[FiJ + 


1 

^6i ^[C][H] 

1-[C]^J^ 


^Z(s) y = s^Z(s) }■ 


(5.23) 
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Typical Section with Zero-Order Compensator 


A zero-order compensator is designed for the typical section of Fig.1 
in incompressible flow. The structural and aerodynamic data sources are 
summarized in Table 2. The minimum-state approximated aerodynamic coef- 
ficients, described in the first minimum state numerical example are used 
to model the system with six structural states and two augmented states. 

The variables of Eq, (5.1) are = [h,a,3,h,a,3,Xai ,Xaz 1 and u = 

3c- The independent variable of the root loci, which defines the flow 
conditions, is the non-dimensional ized velocity V/bwa. The open-loop 
(3c = 0) root loci (Fig. 10) show a violent plunge-torsion flutter at 
V/bwa = 3.02. The control system is required to ensure no flutter be- 
tween v/bwa = 1. and 3.5, and to minimize the cost function of Eq. 

(5.16) as defined in Table 6 for five different design cases. In all 
cases, the process noise is vertical gust velocity Wg, defined by Dry- 
den's PSD of Eq. (3.13) with = 1.0 b^/sec^ and Lg = 50 b. The gust 

forcing vector, ^Ag^ of Eq. (5.13), is derived in Appendix A following 
Bisplinghoff ei [101. 

The closed-loop root loci of cases 3 and 4 of Table 6 are given in 


Fig. 

11. 

The 

mean-square 

response of 3 for a1 1 

the cases is given in 

Fig. 

12. 

The 

mean-square 

« • 

responses of h and a 

for cases 3, 4 and 5 are 


compared to the open-loop response in Fig. 13. 

Assigned complex poles at V/bcJa =3.5 were used as optimization driv- 
ing parameters, one complex pair in the two-measurement cases and two 
pairs in the four-measurement cases. Early analysis, with only one cost 
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TABLE 6 

Typical section control design cases 

Case 

Measurements 

Design 

Cost function 



response 





parameter, Zd 

points, V/bwoi 

weights, ag 

1 

h,h 

3 

2.5, 3.2 

1., 1. 

2 

h, h,a,a 

3 

2.5, 3.2 

1., 1. 

3 

h,h 

3 

2.5, 2.85, 3.2 

1., 4., 1. 

4 

h, h,a,« 

3 

2.5, 2.85, 3.2 

1., 4., 1. 

5 

h,h 

h+3a 

2.5, 2.85, 3.2 

1., 4., 1. 


function point at V/bwa = 2.5, resulted in instability between V/'bwa = 
3.25 and 3.4. This problem was solved by adding V/bwa = 3.2 to the cost 
function points Ceases 1 and 2). As shown in the upper part of Fig. 12, 
the mean-square response of the design parameter 3 at the cost function 
points is reduced by increasing the number of measurements (as ex- 
pected). In between the design points, however, there are response 
peaks, with the four-measurement case being more sensitive to changes in 
flow conditions. 

Adding the response at V/bwa = 2.85 to the cost function (cases 3 and 
4) moves the peaks and the sensitive region to higher velocities 
(lower part of Fig. 12). The root loci of cases 3 and 4 (Fig. 11) show 
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that the two modes which participate in the open-loop flutter mechanism 
(Fig. 10) still interact with each other. But instead of the violent 
flutter, a more complicated interaction, which keeps the system stable, 
is observed. The only instabilities occur at very low velocities as a 
result of the wing-control-surface inertial coupling. This problem is 
easily solved by switching the control system on only after takeoff. 

Figure 13 shows for cases 3 and 4 the mean-square gust responses for 
the translation and pitching velocities, when they are not included in 
the cost function. The closed-loop plunge velocity responses ch^ are 
higher than the open-loop response at subflutter velocities and have 
sharp peaks at higher velocities. The closed-loop pitch rate responses 
are generally lower than the open-loop responses, with the two-meas- 
urement case giving much lower responses than the four-measurement case. 
The main point is that one should keep careful track of response parame- 
ters which are not included in the cost function, 

A change in the design parameter of the two-measurement case to zj = 
h + 3a (case 5) lowers the plunge velocity response oh^ at the high 
airspeed range without changing the subflutter response significantly. 
The improvement in oh^ is accompanied (as expected) by deteriorating 
overall (Fig. 12). 

Research Wing with Zero and Second-Order Compensators 

A compensator is designed to increase the flutter dynamic pressure of 
the research wing of Fig. 4, using a single accelerometer signal and 
minimal control surface mean-square rotation rate. 
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The structural and aerodynamic data sources are summarized in Table 


2. The structural modal properties, as defined in Eq. (2.15), are given 
in Table 7. The control surface is assumed to be structurally uncoupled 
from the wing. Its natural frequency uas arbitrarily chosen to be well 
above the highest uing frequency in this analysis, and the mode shape 
consists of one radian control surface rotation. Open-loop calculations 
by Abel [6] indicate that higher vibration modes do not affect the low- 
est speed flutter mechanism. 


TABLE 7 

Structural modal properties of 

the research 

wing 

Mode 

Natural 
frequency, 
Wn ,i 

(rad/sec) 

General ized 
mass, 

Ms ,ii 

Modal 
Damping, 
? i 

Modal 

deflection 
at sensor 
location, i*! 
(see Fig. 4) 

first bending 

32.88 

1 . kg 

.005 

0.4812 

torsion 

120. 19 

1 . kg 

.005 

-0.2282 

second bending 

161.91 

1 . kg 

.005 

0.2002 

control sui — 
face rotation 

400.0 

.0004 kg-m^ 

.005 

O 

o 


The flight envelope for this design consists of one Mach number, M = 
0.9, and constant air velocity, V = 450 ft/sec (137.16 m/sec), such that 
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the independent variable is the dynamic pressure q. This representation 
fits a wind tunnel test [61. For a given Mach number in atmospheric 
flight* the flight velocity varies someuhat uith q because of the change 
in the speed of sound uith altitude. 

The fourth-order minimum-state approximation of the aerodynamic in- 
fluence coefficients described in the last minimum state numerical example 
is used to construct the plant equation (5.1). The gust coefficients 
^A 3 (ik)^ of Eq. (5.13) are interpolated by Roger's approximation (Eq. 
(4.5)) uith N=6 and y, values of 0.2» 0.4, 0.6 and 0.8. Von Karman's 
pouer spectral density (Eq. (3.14)) represents the vertical gust veloc- 
ity, uith Owg^ = 1.0 m^/sec^ and Lg = 100 ft (30.48 m) . The measurement 
equation (5.2) parameters are defined by Eq. (5.4), uhere the rou vector 
[Hoi takes the values of Table 7. 

The open-loop root loci of the uing modes (Fig. 14) shou an open-loop 
first-bending-torsion flutter at dynamic pressure = 5.37 kPa. Previous 
control designs for this model [6,32] resulted in closed-loop flutter 
dynamic pressures between 9.3 and 9.8 kPa, with compensators of the order 
of 4 or more. The design target here is to achieve similar results with 
a lower order compensator, optimized for minimal cost function Jq = 
at q = 7 kPa. 

The first trial involved the use of a zero-order compensator. Since 
the system is single-input single-output, the only free parameter is a 
single gain C of Eq. (5.7). Closed-loop calculations for various C val- 
ues indicate very limited performance. Tuo typical closed-loop root 
loci are shoun in Fig. 14. A positive control gain reduces both first 
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Figure 14: Open- and closed-loop root loci of the research wing modes 

with a zero-order compensator at M=0.9 
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bending and torsion frequencies with a very small effect on the flutter 
dynamic pressure. A negative control gain decouples the modes of the 
open-loop flutter mechanism but causes torsi on-seoond-bending flutter. 
The variation of flutter dynamic pressure with control gain is shown in 
Fig. 15. Only 105{ increase in flutter dynamic pressure is achieved in 
the very gain-sensitive region of the intersection of the two flutter 
mechanisms, thus this compensator is of little value. 

In the second-order compensator design, the direct control gain is 
set to C = 0. Since Ic = 2 and me = 1» the only free compensator param- 
eters are a^, aj and [Bel (1x2) of Eq. (5.5). The optimal parameters 
are ai = -351.14, az = -28.51 and iBcl* = 1-0.1819,-0.076531. The open- 
and closed-loop root loci of the wing modes are shown in Fig. 16. The 
control system has a considerable effect on the aeroelastic behavior, 
and the flutter dynamic pressure is increased to 9.93 kPa. The compen- 
sator mode interacts with the first bending mode, which gains substan- 
tial damping before fluttering. The torsion and second bending modes 
start to develop a flutter mechanism, as indicated by the turn of the 
second bending mode root locus, but they produce no flutter up to 10 
kPa. 

The variation of op* and op^ of the closed-loop system with dynamic 
pressure is shown in Fig. 17. The response increases monotonical 1 y with 
dynamic pressure, and it does not show the sensitivity of the typical 
section discussed as the preceding example. 
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Figure 15: Variations of the research uing flutter dynamic pressure 

with a zero-order compensator control gain 
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Figure 16: Open- and closed-loop root loci of the research wing modes 

with a second-order compensator at M=0.9 
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PASSIVE FLUTTER SUPPRESSION 


While active flutter suppression systems are still in the theoretical 
and experimental stages, passive flutter suppression means (namely 
structural modification techniques) are well established and widely 
used in the aerospace industry. Optimization procedures to meet flutter 
requirements with minimal weight penalty are given by Markowitz and 
Isakson [40] and by Haftka and Starnes [ 41 ] . These procedures are based 
on repetitive flutter analysis for varying structural parameters. 

An active control system is a method for providing satisfactory 
flutter characteristics with no major structural changes, and thus with 
only a small weight penalty. The designer, however, should not over- 
look local structural modifications as flutter suppression means. Well 
placed concentrated masses or properly tuned springs can suppress 
flutter, as demonstrated in Ref. [42] for wing-store configurations. 

Local passive control means are formulated in the following discussion 
as subcases of active control means, such that a unified active-passive 
analysis can be performed using the methods of the preceding section. 
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Adding a Concentrated Mass 


Although termed as "passive," a concentrated mass mp, added to a 
given structure, is actually a perfect zero-order compensator. It meas- 
ures acceleration and responds with a proportional force. 

The basic assumption of modal analysis is that the motion of the 
structure, before and after the structural changes, is a linear combina- 
tion of the modes which serve as generalized coordinates. This assump- 
tion is never exact for a continuous structure represented by a finite 
number of modes. The analyst should evaluate the assumption's accuracy 
by engineering judgment or by repeating the analysis with a higher num- 
ber of modes. 

In terms of the control terminology used in this paper, the 
measurement is the acceleration at the point where the concentrated mass 
is added. Equation (5.3) becomes 

y = where iHol = ^'i'pf'^ (6.1) 

and where ^$p^^ is the modal deflection row vector at this point. The 
measurement equation (5.2) is constructed using Eq. (5.4). The control 
variable is the reaction force applied by the mass. Equation (5.7) be- 
comes 

u = Cy where C = -mp (6.2) 

The plant equations (5.1) are constructed by leaving the system 
open-loop dynamic matrix [Fi] untouched and using 

(6.3) 
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Uith the additional mass defined as input and output terms of the 
control equations, it can be superimposed on the active control equa- 
tions. The cost function may be supplemented with a weighted mass pen- 
alty term, and the control gain (added mass) can be determined by the 
cost function minimization routine. The additional mass is limited, of 
course, to positive values unless a parasite mass already exists at this 
point. 

The typical section in incompressible flow (see Fig. 1 and Table 2) 
was used to illustrate the effect of adding a concentrated mass. The 
changes in flutter velocity with additional mass at distance Sb real — 
wards of the elastic axis (E.A.) are shown in Fig. 18. Since the elas- 
tic axis is located at a = -0.4, the leading edge is at 8 = -0.6 and a 
reasonable most forward mass location is around 8 = -0.5. The flutter 
velocity can be increased by approximately 20/C with 20/C total mass in- 
crease . 
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Other Passive Control Means 


Other passive means that can be formulated as zero-order compensators 
are springs or dampers connecting two points with modal deflections 
^ and 


For a linear spring with stiffness constant kp» the coefficients of 
the control Eqs. (5.1), (5.2) and (5.7) are 

» iHdl = » iHvl = 0 , 


iHal— 0 » ' IC]— — kp and ICpl — 0 

where ^'Pp^ = ^^pi \ - 


(6.4) 


A viscous damper with damping constant bp is represented by 
^Gz^ = -|»pf » [Hdl = 0 , [Hvl = . 

[Hal — 0 , ^ d — 0 t [C] — —bp and [C^l — 0 


(6.5) 


Unlike the concentrated mass, the spring and the damper have to con- 
nect points with significant relative motion in order to be effective. 

A useful application of pylon stiffness and damping tuning for wing- 
store flutter suppression is given in Ref. [421. 

In some flutter cases, especially when the stability margin changes 
moderately with dynamic pressure, a reasonably-sized dynamic vibration 
absorber can be effective in suppressing flutter. A sketch of a dynamic 
vibration absorber, attached to the wing at a point with displacement Xp 
= is shown in Fig. 19. The equation of motion of mp is 
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C6.6) 


mpZp + bpZp + kpZp = -mpXp 

where Zp = Xp' - Xp 



In state-space form with Zi 


Zp and zz 


Zp, Eq. (6.6) becomes 


^ « 


- 

• 


I" • 


• * 

Z1 


0 

1 


Zi 


0 

, 

• = 





■ + ■ 


Z2 


k p/ni p 

*” b p/m p 


22 


-1 

• ■ 


• 

- 


■ d 




(6.7) 


This is a second-order compensator as defined in Eq. (5.5). The dy- 
namic vibration absorber can be thought of as a mechanical implementa- 
tion of the control compensator and be treated as an active control sys- 
tem. As in the concentrated mass case, the open-loop state vector {z} and 
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dynamic matrix [Fj] remain untouched. The control variable is the 
force applied by the vibration absorber on the structure and the 
measurement is the acceleration at the mounting point. The coef- 
ficient matrices in Eqs. (5.1), (5.3), (5.5) and (5.7) are 


» [Hoi = > 



1 

O 

» 

r ‘’1 



[Ac! = 

It 

0 

CD 


> 

(6.8) 


”kp/iTip "bp/iPpJ 

l-v 




[C] — 0 and [Ccl — Ikp bpl 

The closed-loop equations (5.9) can now be constructed using the matrix 
relations of Eqs. (5.1), (5.4) and (5.8). The design parameters are the 

mounting location, which defines {ipp}, and the absorber parameters m^, b^ 
and kp. The cost-function definition and the optimization procedure can 
be carried out by the methods presented previously for active control, with 
a weighted mass-penalty term added to the cost function. To avoid exces- 
sive relative displacement z^ between the wing and m^, the mean-square 
gust response of can be included in the cost function by adding its 
weighted frequency response to Z(j(im) in Eq. (5.14). The frequency re- 
sponse of Zp is related to the model frequency response of Eq. (5.13) by 

Zp(iw) = •{^p}*'{x(iw) }• (6.9) 

-cj2+i pCJpW+Wp^ 

where Wp = Vk p/mp and ?p - bp/2mpWp 
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The typical section in incompressible flow (see Fig. 1 and Table 2) 
was used to illustrate the effect of a dynamic vibration absorber. The 
changes in flutter velocity with the absorber uncoupled frequency Up, 
for various mounting locations S (see Fig. 18), with mp/ms = 0.2 and Jp 
= 0.2, are shown in Fig. 20. 

When (Op = 0, mp is not connected to the structure and the system is 
unaffected (Fig. 10). When cop •* oo, the flutter velocities approach 
those of the concentrated-mass case (Fig. 18). It is clear from Fig. 20 
that a properly tuned vibration absorber is much more effective than a 
lumped mass. With mp/mg = 0.2, Jp = 0.2 and 8 = -0.5, for example, the 
flutter velocity can be increased by 65%, compared with 20% in the con- 
centrated-mass case. These results, however, are very sensitive to the 
absorber frequency tuning. 
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Wp, rad/sec- 


Figure 20: Variations of flutter speed with vibration-absorber 

frequency, typical section at M = 0. 



CONCLUDING REMARKS 


A review of the basic unsteady aerodynamic derivations and aeroelas- 
tic equations of motion which lead to the widely-used V-g method for 
flutter analysis has been presented. Based on the assumption that the 
wing undergoes simple harmonic motion, the V-g method gives correct re- 
sults only for the flutter boundary flow conditions. Matched point 
solutions must be determined by Iteration with density as a parameter. 

Direct stability margin calculations at pre- and post-flutter condi- 
tions need aerodynamic Influence coefficients for arbitrary motion. 

In order to solve the stability equations by the methods of linear 
algebra, the influence coefficients are approximated by rational func- 
tions of the nondlmensionalized Laplace variable s'. Such approxima- 
tions lead to constant coefficient, linear, finite state-space aeroelas- 
tic models. 

Several rational approximation techniques are discussed, and a new 
"minimum-state" approximation method, which yields a minimal-order aero- 
elastic model for a given accuracy, is presented. All the methods con- 
sist of aerodynamic data calculated for oscillatory motion. The mini- 
mum-state method is first applied to a three-DOF typical section in in- 
compressible flow. With two aerodynamic augmented states, the 
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approximated influence coefficients are within 5% of the "exact" 
Theodorsen coefficients in a range of i30° of the imaginary axis of the 
s' plane. Root loci of important structural modes usually fall within 
this range. 

Application of the minimum-state method to the typical section at M = 
0.7 and to the four-mode research wing model shows that, when the number 
of aerodynamic augmented states is similar to the number of degrees of 
freedom, all the approximated coefficients are within 5% of the exact 
ones along the imaginary axis. Off the imaginary axis there were no 
available data with which to compare the results . The minimum- 
state approximation is demonstrated to be significantly superior to 
Roger's and the matrix Fade approximants in terms of accuracy per model 
order. One disadvantage of the minimum-state method is that it is more 
complicated and computer time-consuming, because it involves solving 
nonlinear least-squares problems. This investment, however, pays off in 
repetitive flutter and control calculations. 

It is shown that the minimum-state aeroelastic model can be used to 
design constant-parameter active control systems which by minimizing a 
cost function, defined as mean-square gust response over several flight 
conditions, assure stability over the entire flight envelope. Numerical 
examples for a typical section with a zero-order compensator and a 
research wing with a second-order compensator demonstrated the design 
procedure and showed a significant favorable change in the aeroelastic 
behavior with a single trailing edge control surface. 
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Finally, more traditional passive control means — such as concentrated 
masses, springs, dampers and dynamic vibration absorbers — are formulated 
as control input and output terms. They are shown to be analogous to ac- 
tive control means, making it possible to carry out simultaneously 
active and passive control analyses. 


RECOMMENDATIONS 

1. The minimum-state approximation should be applied to more 
aeroelastic systems in order to explore further its merits as 
compared with other rational approximation methods. The com- 
putation time might be reduced by finding more efficient solu- 
tions to the nonlinear least-squares problem stated in Eq. 
(4.27). Other improvements may be obtained by allowing the 
aerodynamic "open-loop" roots (the eigenvalues of [R'D to be 
complex , 

2. Unsteady aerodynamic theory and computational algorithms 

should be further developed to provide aerodynamic influence 
coefficients in the entire s' plane for three-dimensional 
wings in compressible flow. This will enable checking the ac- 
curacy of the rational approximations off the imaginary axis 
of the Laplace domain. The minimum-state procedure should be 
modified to accommodate data for arbitrary s' values. 

3. Rational approximations of the frequency-dependent gust loads 
and "coloring" of the gust velocity power spectral density 
should be attempted. This will add augmented states to the 
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aeroelastic model but uill also enable replacement of the 


heavily time-consuming numerical integrations for the mean- 
square gust response by algebraic solution of Lyapunov equa- 
tions 1491. 

4. The suggested active control design procedure should be com- 
pared with other fluttei — suppression control techniques by ap- 
plying them to the same mathematical models, wind-tunnel tests 
and flight test demonstrations. 

5. The applicability of the control models developed in this study 
are not limited to the design of constant parameter, low-order 
continuous control systems. These models can be used in the 
design of other control systems, such as scheduled parameter 
systems, and in controllability, observability and sensitivity 
studies for choosing the control devices. The relatively 
low-order model has special potential for developing real- 
time, adaptive, digital control systems in which computation 
time is a critical parameter. 

6. The merits of formulating passive flutter-suppression means as 
input and output control terms should be further explored by 
combining them with active control design and by applying them 
to the design of special flutter-related devices, such as de- 
coupler pylons, control -surf ace mass balance and control-sys- 
tem actuators. 
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APPENDIX A 


UNSTEADY AERODYNAMICS OF A TYPICAL SECTION IN INCOMPRESSIBLE FLOW 

This appendix summarizes the unsteady aerodynamic loads on a typical 
section in incompressible flou due to arbitrary section motion and due 
to a disturbing gust with sinusoidal vertical velocity. 

The unsteady aerodynamic loads on the typical section of Fig. 1 in 
tuo-dimensionl incompressible flou due to oscillatory section motion 
were derived by Theodorsen [12] and generalized to arbitrary motion by 
Edwards [19], following Sears [131. The Laplace transform of the loads 
reads 

^LaCs)[ = q[A(s")]]X(s)^ (A.1) 

where •{xCt)^'^ = [h/b,a,B] , ■{LaCt)}-* = [-Lb,Ma»MoJ, 

q is the dynamic pressure, s' is the nondimensional Laplace variable 
sb/V, and [ACs'J] is the generalized aerodynamic influence coefficient 
matrix 


[A(s')l = 2b2 


(Mncls'^ + 


[Bncl+C(s')^Rl [[Szl S' 


+ [Kncl+C(s')^Ri [[Si ] 


(A. 2) 


where the coefficient matrices are 




—If 

ira 

Ti 

iMncl - 

ira 

-iT(a^+ 1 / 8 ) 

“Ti 3 


Ti 

“Ti 3 

T3/TI 


0 

-IT 

Th 

[BncI = 

0 

Tf(a-i) 

“Tl 6 


0 

-Ti 7 

-Ti 9/11 


0 

0 

0 

[ K PC ] - 

0 

0 

-Ti 5 


0 

0 

-T 1 Q/"a 

= I- 2 tt 

2ii(a+5) 

-Ti 2 1 


[Si ] = [ 0 1 1 

IS 2 I = 11 l-a ] 

and the constants are 

Ti = -(1/3) (2+c2)/r^ + c cos-’c 

Ta = -(1/8)(1 -c2)(5c 2+4) + C 1/4)c ( 7+20^ )i/ l-c^ cos'^c 
- (c^+ 1/8) (cos' ’c)2 

Th = ct/I-c^ = cos'*c 

Tj = -d-c^) - (cos*’c)^ + 2 ct/i-c^ cos’’c 
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T? “ ( 1/8) c ( 7+2c^)i/l"C^ - (c^+1/8) cos'*c 
Ts = -(1/3)(1+2 c2)a/i-c 2 + c cos-’c 
T, = H ( V3)(1-c2)3''2 + aT^] 

Tio = t/T^ + cos’’c 

Tii = (2-c)\/l-c2 + (1-2c) cos'^c 

Ti 2 = t/T-^(2+c) - (2c+1)cos‘’c 

Ti 3 = -ilT; + (c-a)TiI 

Ti 5 = T^ + Tio 

Tie = Ti - Te - (c-a)Ti| + iTn 
Ti7 = -2T, - Ti + (a-i)T4 
Tie - Ts - T^T1o 

» 

Ti9 = -^T^T11 

The generalized Theodorsen function is 
Ki (s') 

C(s') = (A. 3) 

Ko(s')+Ki (s') 

where Kn is a modified Bessel function of order n. Using Eq, (9.6.4) of 
Ref. [ 50 J for 0 < arg(s') < tv, Eq. (A. 3) becomes 

Hi'2>(z) 

C(s') = where z = -is' (A. 4) 

Hi < 2 ^ (z) + i Ho^ ^ ^ (z) 
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where is the Hankei function which is re?ated to Bessel functions 

by 


Hn<2’(z) = Jp(z) - iYp(z) 


(A. 5) 


As indicated by Edwards [19l> as |s'| 0, C(s') -» 1 independent of 

approach direction. To calculate C(s')» ascending power series of Bes- 
sel functions are derived from Eqs. (9.1.10) to (9.1.13) of Ref. [501. 


j 

CO (-z2/4) 

Jo(z) = 1+5; (A. 6) 

j=1 (j!)2 




2 

' 

2r 00 r(-22/4) j -l' 

Yo(z) = - 

jHn(z/2)+v 

Jo(z) - - Z ( Z 1/m) 

TT 

■ ^ 

TT>.j = lL (j!)Z m=i -1. 


where y = 0.577215665 


Jl (z) 


2 ' 
- 1 
2 - 


3 

CO (-2^/4) I 

- X 


j=1 j! ( j+1)!J 


(A. 8) 


Yi (z) 


2 ■ 
IT . 


- + lJln(z/2)+ylJi (z) 
z 


z ' 
— 1 
2ti - 


00 r j 1 1 (-2^/4) ' 

+ Z 2(^ Vm) + 

j = 1 *- m=1 j+P j ! ( j+1 ) !- 


(A. 9) 


The ascending series are converging very fast for the study range of in- 
terest 0 < |s'| 1 2. 
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The unsteady aerodynamic loads due to a disturbing gust uith 
sinusoidal vertical velocity of amplitude Ug" and frequency a read 

^Lg(ik,t)[ = (q/V)^ Ag( ik ) ^Ug"exp C iwt ) (A. 10) 


where k= wb/V. The hinge and pitch components of the gust loads are 
given by Eq. (5-375) of Ref. 110]. In our notations, the related compo- 
nents of •{Ag^ read 


Ag.i = -47Tb2^C(ik)(Jo(k)-iJi(k)] + iJi(k)[ 
Ag ,2 - ~ ^ J + a)Ag ,1 


(A. 1 1) 


For a relatively stiff actuator, and when the open-loop gust-induced 
hinge moments are not of particular interest, their effect is negligible 
and it is assumed that Ag ,3 = 0 . 
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